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Abstract 


The formation of structures in living organisms is addressed within the framework of far- 
from-equilibrium chemical systems using methods of statistical physics, such as kinetic 
theory and stochastic methods, at an intermediate, mesocopic scale. Three directions are 
explored. For the purpose of investigating the stochastic elimination of a fast variable, a 
fast species is eliminated from a nonlinear chemical mechanism. The fluctuations of the 
slow species using Langevin equations and a master equation are not correctly predicted 
by the reduced mechanism. The coupling between the fluctuations and the nonlinearities 
of deterministic dynamics makes the use of the quasi-steady-state approximation deli- 
cate when the studied system requires a good control such as in fluorescence correlation 
spectroscopy (FCS). A submicrometric Turing pattern is simulated in a concentrated sys- 
tem in order to refute certain objections to Turing’s model regarding the preservation 
of proportions in embryos. Assuming an appropriate role of the solvent in the chemical 
mechanism is sufficient to control the wavelength of the structure by monitoring the con- 
centration of the solution. The results can be exploited to design materials with controlled 
submicrometric properties in chemical engineering. Following a biomimetic approach, ex- 
perimental conditions leading to the termination of the Turing structure associated with 
a decrease of the wavelength are proposed. The sensitivity of the Fisher-Kolmogorov, 
Petrovsky, Piskunov wave front to small perturbations is used to characterize the effects 
of the deviation from the dilution limit on diffusion. As a result, the shift of the concentra- 
tion profiles of two species associated with different diffusion coefficients is a well-adapted 
criterion to detect perturbations induced by high concentrations. Contrary to the results 
of a deterministic description, the front speed deduced from the master equation in the 
dilute case sensitively depends on the diffusion coefficient of the consumed species. In the 
case of a concentrated solution, the properties of the wave front obtained in the dilute 
case remain valid but are mitigated by cross-diffusion terms which reduce the impact of 
different diffusion coefficients. 
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Streszczenie 


Tworzenie struktur w organizmach żywych rozważane jest w ramach odległych od równowagi 
układów chemicznych przy użyciu metod fizyki statystycznej, takich jak teoria kinetyczna 
i metody stochastyczne, w pośredniej, mezoskopowej skali. Badane są trzy kierunki. W 
celu zbadania eliminacji szybkiej zmiennej stochastycznej, wprowadzono szybko reagu- 
jący związek do nieliniowego mechanizmu chemicznego. Fluktuacje związku o powolnej 
dynamice uzyskane za pomocą równań Langevina i równania master nie są prawidłowo 
przewidywane w mechanizmie zredukowanym. Sprzężenie fluktuacji z nieliniowością dy- 
namiki deterministycznej sprawia, że stosowanie przyblizenia quasi-stacjonarnego jest de- 
likatne, gdy badany układ wymaga dobrej kontroli, np. w spektroskopii korelacji fluo- 
rescencyjnej (FCS). Submikrometryczna struktura Turinga jest symulowana w układzie 
stężonym w celu odrzucenia pewnych zastrzeżeń do modelu Turinga dotyczących zachowa- 
nia proporcji w zarodkach. Przyjęcie odpowiedniej roli rozpuszczalnika w mechanizmie 
chemicznym jest wystarczające do kontroli długości fali struktury poprzez monitorowanie 
stężenia roztworu. Wyniki mogą być wykorzystane do projektowania materiałów o kon- 
trolowanych właściwościach submikrometrycznych w inżynierii chemicznej. Zgodnie z 
podejściem biomimetycznym, proponowane są warunki doświadczalne prowadzące do za- 
kończenia struktury Turinga związane ze zmniejszeniem długości fali. Wrażliwość frontu 
fali Fishera-Kolmogorowa, Petrovskiego, Piskunova na małe perturbacje jest wykorzysty- 
wana do scharakteryzowania wpływu odchylenia od granicy rozcieńczenia na dyfuzję. W 
rezultacie, rozsunięcie profili stężeń dwóch składników związanych z różnymi współczyn- 
nikami dyfuzji jest kryterium dobrze dostosowanym do wykrywania perturbacji wywołanych 
przez wysokie stężenia. W przeciwieństwie do wyników opisu deterministycznego, pręd- 
kość frontu wyprowadzona z równania master w przypadku rozcieńczonym zależy od 
współczynnika dyfuzji gatunku konsumowanego. W przypadku roztworu stężonego, właś- 
ciwości frontu fali uzyskane dla przypadku rozcieńczonego pozostają ważne, ale są łagod- 
zone przez efekty dyfuzji krzyżowej, które zmniejszają wpływ odmiennych współczyn- 
ników dyfuzji. 
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Preamble 


Living organisms are fascinating examples of reaction-diffusion systems evolving into self- 
organized structures through sustained exchanges with the environment. Periodic time 
oscillations of concentrations, chaotic behaviors, and complex spatial structures are ob- 
served in chemical systems maintained far from equilibrium [1, 2, 3, 4]. 

The formation of Turing patterns and the propagation of a wave front are recognized 
as playing an essential role in the structuring of living organisms [2]. In 1952, Tur- 
ing proved that reaction-diffusion processes in initially homogeneous far-from-equilibrium 
systems may lead to the formation of periodic spatial patterns observed in biology [5]. 
Zebra stripes and leopard spots are classical illustrations of the contribution of far-from- 
equilibrium reaction-diffusion systems to the modeling of morphogenesis. Turing’s model 
is now used as a chemical basis for embryo development, e.g. in limb formation and teeth 
development [6, 7]. Even earlier, in 1937, Fisher built a model describing the propagation 
of a favored trait in a population. Simultaneously, Kolmogorov, Petrosky, and Piskunov 
studied the traveling solution of the same equation, further referred to as the FKPP wave 
front. The information conveyed by a signaling wave front during segmentation is at the 
base of different models of development [8, 9, 10, 11, 12]. 

Although associated with chemical mechanisms involving elementary steps between 
molecules, Turing patterns and FKPP wave fronts were first studied within the framework 
of a macroscopic approach based on partial differential equations. However, the descrip- 
tion of phenomena arising during the early development, when the embryo is composed of 
a small number of cells, may require approaches at smaller scales. The interplay between 
fluctuations and nonlinear dynamics is known to induce non intuitive, model-specific be- 
haviors [1, 13, 14]. The description of reaction-diffusion systems at the mesoscopic scale 
requires stochastic methods introducing random variables but still ignoring the detail of 
molecular dynamics [15, 16]. The crudest stochastic method used to describe reaction- 
diffusion systems consists in adding a Langevin force to the deterministic equations and 
assuming that the probability distribution of the concentrations is Gaussian [13, 14]. 
The correct description at the mesoscopic scale leads to a master equation relying on 


the well-founded hypothesis that reactions and diffusion are Markov processes. Reac- 
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tions are seen as birth and death processes, whereas diffusion is interpreted as jumps 
between adjacent spatial boxes [1]. Numerical simulations become necessary to describe 
dynamics at the molecular scale. When the details of electronic reorganization within a 
molecule can be ignored, classical mechanics is sufficient to describe particle dynamics. 
Simulations of molecular dynamics involve the computationally expensive deterministic 
processing of particle displacements and collisions, which could compromise reaching the 
space and time scales necessary for the emergence of structures in growing systems. In 
these conditions, it is appealing to consider the smart method intuited and developed 
by Graeme Bird and known as the direct simulation Monte Carlo (DSMC) [17, 18, 19]. 
DSMC has been designed to simulate the dynamics of dilute gases and is particularly 
adapted to space applications [20, 21, 22]. The collisions are efficiently treated through 
a random procedure known as Monte Carlo. It has been proven [23] that DSMC simu- 
lates the Boltzmann equations governing the evolution of the distribution functions for 
position and velocity of the particles [24]. In addition to following the laws of kinetic 
theory, DSMC provides stochastic trajectories and correctly simulates the fluctuations 
in agreement with the master equation [25]. Hence, the direct simulation Monte Carlo 
method offers an efficient alternative to molecular dynamics and gives access to space and 
time scales compatible with the simulation of emerging micrometric structures. 

In this work, I developed stochastic approaches to far-from-equilibrium structures and 
focused on Turing patterns and Fisher-Kolmogorov, Petrovsky, and Piskunov wave fronts, 
both for their relevance in biology and the richness of their behavior from a fundamental 
point of view. Within this framework, several usual approximations have been revis- 
ited. The question of the elimination of a variable, encountered in the study of Turing 
patterns, incited me to investigate the validity of the widely used steady-state approxi- 
mation in systems with large fluctuations. I was also led to consider the deviation from 
the high-dilution limit as a possible way to tune the features of a pattern. Beyond the ap- 
plication to the adaptability of a structure, dealing with concentrated systems prompted 
me to deepen my knowledge of the cross-diffusion phenomenon and the associated form 
of Fick’s law deduced from irreversible thermodynamics in the linear domain [26]. In par- 
allel, extending both the master equation and the direct simulation Monte Carlo method 


to concentrated solutions was an attractive challenge. 


The manuscript is organized as follows. In the Chapter I, I recall the analytical and 
numerical methods that I used and adapted to concentrated systems. 

The validity of the steady-state approximation in a small chemical system is questioned 
in Chapter II. In order to investigate to which extent the description of the fluctuations 
remains correct after elimination of a fast variable, I compared the correlations of concen- 


tration fluctuations for two different chemical mechanisms leading to the same reduced 
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mechanism in which a reactive intermediate species has been omitted. I first drew an 
analogy between the steady-state approximation and the Born-Oppenheimer approxima- 
tion and expressed the conditions of validity of the steady-state approximation applied at 
the macroscopic scale. Then, I developed two stochastic approaches. Using an analytical 
Langevin approach and simulations of the master equation, I showed that the correlations 
of concentration fluctuations of the slow species are different depending on whether a fast 
intermediate species is considered or not in the reaction scheme. In biology, Fluores- 
cence Correlation Spectroscopy (FCS) is widely used to study the dynamics of labeled 
species, for example to evaluate rate constants when the reaction scheme is supposed 
to be known [27, 28]. The interpretation of the results requires the comparison of the 
experimental data with analytical expressions of the correlations. My results point out 
that, even if it enables the analytical computation of the correlations, a tractable reduced 
reaction scheme could be misleading. The results have been published in G. Morgado, 
B. Nowakowski, and A. Lemarchand, Elimination of fast variables in stochastic nonlinear 
kinetics, Phys. Chem. Chem. Phys. 22, 20801 (2020) [29]. 


Chapter III is devoted to Turing patterns. After recalling the basics of Turing struc- 
tures, I place the subject in the context of morphogenesis. Turing’s model relies on a 
remarkably small number of processes involving two initially homogeneously distributed 
chemical substances that interact to produce stable patterns. The model involves two 
chemical species, an activator and an inhibitor. The minimal reaction-diffusion scheme 
for the emergence of Turing patterns requires the autocatalytic production of the acti- 
vator and the faster diffusion of the inhibitor [30]: The structure develops through local 
self-enhancement in conjunction with long-range lateral inhibition [31]. The wavelength 
of the periodic spatial structure is determined by the reaction rate constants and the 
diffusion coefficients of the chemical species. Contrary to Taylor vorticies in hydrody- 
namics, the striking property of Turing patterns is that the wavelength of the structure 
does not depend on the boundary conditions. The robustness of Turing patterns is a 
strong feature, but also an argument against them in morphogenesis. Indeed, Turing pat- 
terns lack scaling properties: they do not account for size adaptation of the wavelength 
to the global size of the system. Yet, models of somitogenesis should reflect that the size 
of the vertebrae is proportional to the size of the embryo [32, 33, 34, 35]. I addressed two 
points related to the growth of a Turing pattern in a growing system, the question of the 
scaling properties of a periodic spatial structure and the question of the termination of 


the structure. 


Recently, the Polish-French group proposed to solve the problem of scaling of a Turing 
pattern at the macroscopic scale by considering the possible perturbations induced by high 


concentrations of reactants [36]. In these conditions, the variation of the concentration 
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of the solvent cannot be neglected. The model that the group developed includes the 
participation of the solvent into the reaction scheme. This third substance introduces an 
additional variable concentration that can be harnessed to control the behavior of the 
system. Changing the dilution of the medium makes it possible to tune the wavelength 
of the emerging Turing pattern. During this PhD, I performed a deeper analysis of 
the three-variable model at the microscopic scale and examined if the deviation from 
the high-dilution limit also induces a control of the pattern in small systems. A proof 
of concept based on simulations of particle dynamics was necessary. To this goal, I con- 
sidered the direct simulation Monte Carlo (DSMC) method and adapted it to concentrated 
solutions. The simulations show that doubling the concentration of the solute leads to 
decreasing the wavelength of the structure by a factor of 2. The results can be considered 
as a possible interpretation for proportion preservation of embryos in morphogenesis. 
They can also be used to design materials with controlled submicrometric properties in 
chemical engineering. The results have been published in G. Morgado, B. Nowakowski, 
and A. Lemarchand, Scaling of submicrometric Turing patterns in concentrated growing 
systems, Phys. Rev. E 98, 032213 (2018) [37] and were presented at the 31st International 
Symposium on Rarefied Gas Dynamics in Glasgow, in 2018 [38]. 


The question of the termination of a spatial structure is compelling in morphogenesis: 
The spine of the vertebrates ends with smaller vertebrae and the fingers with smaller 
phalanges. Within the framework of Turing patterns, this phenomenon implies both a 
decrease of the amplitude of the spatial oscillations and a decrease of the wavelength. 
Deciphering the mechanisms actually controlling the termination of the spine in an em- 
bryo is far beyond the scope of this work. The aim was to propose a possible mechanism, 
inspired by biological structures and compatible with the implementation in a chemical 
engineering context. The boundary conditions chosen by the group in 2016 [36] to model 
the growth of the spine behind a freely propagating wave front are well adapted to the 
design of an artificial spatial structure. I performed a systematic analysis of the effect 
of all rate constants and diffusion coefficients on the stability and the wavelength of the 
structure. Interestingly, a monotonous variation of almost any of the dynamical param- 
eters leads to the simultaneous loss of stability of the structure and the decrease of the 
wavelength. Only the variation of the diffusion coefficient of the activator causes anticor- 
related results. Locally varying a rate constant or the diffusion coefficient of the inhibitor 
in a given chemical system is not straightforward from an experimental point of view. 
For an easy implementation in chemical engineering, I suggest to impose an appropriate 
spatial profile for the concentration of the reservoir of inhibitor, resulting in the expected 
variation of the effective rate constant controlling the injection of the inhibitor in the 
system and leading to the desired termination of the structure. The results have been 


published in G. Morgado, L. Signon, B. Nowakowski, and A. Lemarchand, Termination 
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mechanism of Turing patterns in growing systems, Acta Phys. Pol. B, 50, 1369 (2019) [39]. 


The results obtained for the model of Turing patterns with a reactive solvent have 
opened new directions of research. The deviation from the high-dilution limit has an 
impact on both reaction and diffusion. In the case of Turing patterns, the perturbation of 
diffusion induced by high concentrations has a negligible effect on the wavelength of the 
structure. In order to investigate the possible consequences of the perturbation of diffusion 
in a concentrated reaction-diffusion system, I considered the case of pulled wave fronts, 
known to be sensitive to even small perturbations [40]. The results dealing with wave 
fronts are given in Chapter IV. For species with identical diffusion coefficients, the group 
already showed that an FKPP front is sensitive to the discrete nature of particles [41] 
and to reaction-induced deviations from partial equilibrium [42]. I proposed an FKPP- 
based model involving two species A and B engaged in the reaction A + B —> 2A 
with different diffusion coeffcients. In a concentrated system, the resulting wave front 
turns out to be a sensor revealing perturbations of diffusion at the macroscopic scale. 
Specifically, I showed that the difference of concentrations between the two species A and 
B at the inflection point of the A profile is a good indicator for diffusion perturbation in 
concentrated systems. The results have been published in G. Morgado, B. Nowakowski, 
and A. Lemarchand, Fisher-Kolmogorov-Petrovsky-Piskunov wave front as a sensor of 
perturbed diffusion in concentrated systems, Phys. Rev. E 99, 022205 (2019) [43]. 

Deviations from the high-dilution limit are more prone to happen in small systems, 
typically in a living cell, where the amplitude of concentration fluctuations are signifi- 
cant. I therefore performed a stochastic analysis of fluctuation effects on an FKPP front 
with perturbed diffusion in a mesoscopic system. Unexpected results on a more than 
eighty-year-old problem have been obtained: In a dilute system of small size, the wave 
front propagates more slowly than expected if species B diffuses faster than species A. 
In a concentrated system, the results are mitigated by cross-diffusion which reduces the 
effect of different diffusion coefficients. The results have been published in G. Morgado, 
B. Nowakowski, and A. Lemarchand, Stochastic approach to Fisher and Kolmogorov, 
Petrovskii, and Piskunov wave fronts for species with different diffusivities in dilute and 
concentrated solutions, Physica A 558, 124954 (2020) [44]. 


Chapter V contains conclusions. 
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Chapter I 


Methods 


This chapter presents different theroretical tools of macroscopic and stochastic descrip- 
tions of dynamical systems in chemistry that will be used in later chapters. Mathematical 


notations are also introduced. 


1.1 Chemical kinetics 


All considered systems involve reactive species and possibly non reactive species such as 
the solvent. The steps of a reaction scheme are supposed to be elementary reactions for 
which the standard laws of kinetics apply. The left-hand side of elementary steps involve 
one or more molecules. When only one molecule is involved in the left-hand side, the step 
leads to a first-order reaction rate. At the microscopic scale, it usually corresponds to 
an internal molecular reorganization. When two molecules are involved in the left-hand 
side, the step leads to a second-order reaction rate. At particle scale, it occurs through 


reactive collisions between two molecules. 


It is very unlikely that a collision of more than two molecules occurs. Therefore, 
we consider that reactions of order higher than two result from the reduction of a set 
of elementary steps of first and second orders. The assumptions making this reduction 
possible will be discussed in Sec. II. Each chemical species is assumed to be subject to 


diffusive transport. 


For some systems, we introduce reservoirs. A reservoir maintains the concentration 
of a chemical substance constant by instantaneously removing or adding molecules when 
needed. A reservoir is denoted by the letter R in a reaction scheme, and its constant 


concentration denoted by cr. 
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2 Chapter I. Methods 


1.1.1 Rate equations 


We consider a reaction scheme involving m different steps with n different species X;. 


Each step 7 is associated with a rate constant k; such that 


DX — Y Bij; (1.1) 
i=l i=l 


where a; ; and f; ; are possibly vanishing stoichiometric coefficients. 
The rate equations for the concentrations c; of species X; associated with this reaction 


scheme are 


dic; = 5 kj [Bij = Qi j] I I ee (1.2) 
j=l 


where the symbol d; denotes the derivative with respect to time ap 
The system state at a given time ¢ is then defined by the vector of concentrations 
c = (c(t), c2(t), ..., Cn (t)) 
If the chemical species X; is involved in a second or higher-order reaction, the system 


given in Eq. (1.2) is nonlinear and, in most cases, has no analytical solution. 


1.1.2 Reaction-diffusion equation and Fick's law 


In addition to the reaction, we introduce diffusive transport. Diffusion tends to homoge- 
nize concentrations in an inhomogeneous medium. According to Fick’s law, the diffusive 


flux ji of a given chemical species X; is proportional to the concentration gradient 
ji = —DiVci (1.3) 


where D; is the diffusion coefficient of species X,. In a reaction-diffusion system, the local 
variation of concentration c; of species X, is expressed as the sum of a reactive term and 


a diffusive term, where the latter is given by the divergence of the flux 


n 


dici = Y kj (Big — ois] [] = Ve (1.4) 
j=l 


il-1 
A kl 0 
where the symbol ô; denotes the partial derivative with respect to time ar 


In concentrated mixtures involving more than two species, cross-diffusion effects may 


appear and are discussed in Sec. 1.3 
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Equation (1.4) is valid in a macroscopic approach, in which local fluctuations have 


been neglected. 


1.1.3 Linear stability analysis 


In the general case of a homogeneous chemical system, dynamics is governed by a nonlinear 


system of differential equations for the concentration vector c = (c1, ..., Cn) 
dc = f(c) (1.5) 


where f is a vector function characterizing the reaction rates. 


The local dynamics around a steady state c' = (cf, ..., 2) obeying 


dic? = 0, Vi (1.6) 


can be studied by the linearized dynamics around this state. In order to perform a linear 
stability analysis, we introduce the small deviation dc; = c; — c? to the steady state which 


obeys 
dzóc = Mc (1.7) 


where dc = (dcq,...; Öcn) and 


is the Jacobian matrix for c = c", called the stability matriz. Some laws of conservation 
can be observed and reduce the number of independent variables. We consider that the 
n variables are all independent. The deviation dc is related to the vector y = (41, -.., Yn) 


in the eigenbasis of M? by 
ôc = Py (1.9) 


where P is the change-of-basis matrix. The linear differential equations given in Eq. (1.7) 


lead to uncoupled equations in the eigenbasis of MY 


diyi = Hili (1.10) 


where u; are the eigenvalues of MY. Equation (1.10) is straightforwardly solved, leading 
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Fig. I.1 The five typical phase portraits for an n = 2 system. [2] 
to 
yi = y erit (L11) 
where 44"! is the initial value of y. 


From Eq. (1.11), it appears that the linear stability of the system is governed by the 
eigenvalues u;. If the real parts of all eigenvalues are negative, then the system eventually 


0. The time 7; = 1/|Re(pui)| corresponds to the 


converges towards the steady state c 
typical relaxation time in the direction associated with 7;, as long as the deviation dc is 


in the linear domain around the steady state c'. 


In Fig. 1.1, different phase portraits are presented for an n = 2 system. Five typical 
phase portraits exist: the stable node corresponds to two real negative eigenvalues, u1 < 0 
and u2 < 0, the stable focus corresponds to two complex-conjugate eigenvalues whose real 
parts are negative, Re(u1)£ 0 and Re(øu2)< 0, the unstable node corresponds to two real 
positive eigenvalues, u1 > 0 and ua > 0, the unstable focus corresponds to two complex- 
conjugate eigenvalues whose real parts are positive, Re(i4)> 0 and Re(u2)P 0, and the 


saddle corresponds to two real eigenvalues with different signs, det M? = pu; a < 0. 
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Finally, substituting Eq. (1.11) for % into Eq. (1.9), we obtain 


de = 5 Par erit (1.12) 
j 


Linear analysis describes local evolution around a steady state within the framework 
of a macroscopic, deterministic approach. In small systems, close to situations where the 
dynamics is sensitive to small perturbations, such as bifurcations, a deterministic analysis 
may be insufficient. A stochastic approach, including the description of the fluctuations 
at the mesoscopic scale, is then required. 

I have used linear stability analyses extended to inhomogeneous systems to study the 


termination of a Turing structure [39] presented in Sec. III.2. 


1.2 Stochastic chemical kinetics 


Although chemical dynamics is driven by discrete elementary processes, it is usually 
sufficient to consider deterministic equations to describe the macroscopic evolution of 
a chemical system. However, fluctuations may not be negligible in small systems and a 
stochastic approach may be required [1, 16, 45]. In this section, we introduce two different 
stochastic descriptions of a chemical system, the chemical Langevin equations [13, 14] and 


the master equation [1]. 


1.2.1 Chemical Langevin equations 


In this subsection, we introduce a probabilistic approach to a reactive system [46]. We 
consider the vector of concentrations c as a vector of random variables. Intuitively, the 
propensity or transition rate p;(c) that the j-th step of the reaction occurs in the next 
time interval [t, t + dt] depends on the order of the step. The probability for a first-order 
step to occur is proportional to the number of molecules N; in the system, since each 


molecule is susceptible to be re-organized 


pi (c)dt =Probability that a molecule re-organizes itself 


x Number of molecules N; (1.13) 
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The probability for a second-order step to occur depends on the product of two probabil- 


ities 


pi (c)dt =Probability that two given molecules collide 
x Probability that two colliding molecules react 


x Number of possible pairs of molecules (1.14) 


According to kinetic theory, the first probability is proportional to the average relative 
speed and collision cross-section of the two molecules and inversely proportional to the 
system size (2. The second probability expresses that reaction occurs if the collision 
energy exceeds a certain threshold known as activation energy. The product of these two 
probabilities gives the so-called rate constant introduced in Sec. I.1, except that specific 
care is needed in order to convert a probability derived from a discrete number of molecules 
into a proportionality factor that deals with continuous concentrations. Specifically, the 


conversion introduces as much (2 factors as the order of the reaction. 


If concentrations c; are locally homogeneous, p;(c) is typically similar to the reactive 


term in Eq. (1.2) for all reaction orders 


n 
parle) = OG, I I gå (1.15) 
i=1 
During a finite but small time interval [t,t + 7], the number of reactions r; (c,7) of step 


j is a random variable whose mean (r;);,, is deduced from Eq. (1.14) 


(Tz)t,r = pj(e)T (1.16) 


where the concentration c is evaluated at time t. For the reaction scheme given in Eq. (1.1), 


the concentration c; at time t + 7 is given by 
1 m 
cilt + 7) = c(t) + 52, (Big — 04,5) ri (6,7) (1.17) 
j=l 


Some conditions must be fulfilled for Eqs. (1.16) and (1.17) to hold. On the one hand, 7 
must be sufficiently small for the variations of concentration between two consecutive time 
steps to be small. It implies that the propensity given in Eq. (1.14) is constant over the 
time interval [t, t + 7]. This condition is satisfied if the number of each type of molecules 
in the system is much larger than 1. On the other hand, 7 must be sufficiently large for 
rj (c, T) to be substantial, i.e. for the mean number of reactions (r;)+, to be much larger 
than one. It is not unusual to find systems with sufficiently large numbers of molecules 


that respect both conditions. Typically, mesoscopic systems are of adequate size for these 
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two conditions to be satisfied. 


It can be argued that each random variable r; (€,7) follows an independent Poisson 
distribution of mean u. However, the construction of a standard Langevin equation 
introduces independent first and second cumulants of the probability distribution of c. 
Using the condition that the system contains a large amount of molecules, each Poisson 
random variable r; (c,7) can be approximated by a normal random variable N (p, o?) of 


same mean u and variance o”. The linear combination theorem for normal distributions 
N (u, 07) = p+ 0 N(0, 1) (1.18) 
and Eq. (1.16) allow us to write Eq. (1.17) into the form 


1 m 
alt +T) = c(t) + TÈ (Bij — Qij) [p;(c)r ES (p;(e)r)" PN (0, 1)] (1.19) 
2 
Considering the time 7 as an infinitesimal time interval dt that respects the conditions 


mentioned above and using Eq. (1.15), we write Eq. (1.19) as a chemical Langevin equation 


m n n 1/2 
> AS Qi! Qi! 
= pa (Bij — Qij) He y + aj Bij — Qij) E C. i å Ej (t) (1.20) 
j=l =! U=l 
with independent Gaussian white noises 6; (ż) 
(5; (Ł)) = 0 
(E (DES) = djal — 1”) (1.21) 


In Eq. (1.20), the first term is the deterministic rate equation given in Eq. (1.2) and 
the second term is a noise term denoted by mi. The noise 7; is written as the sum of the 


noises 1); ; associated with the reaction steps j involving the chemical species X; 


m m n 1/2 
= 20) > (bij — Qij) E ġo) ċj(t) (1.22) 


U=l 


The variances and covariances of the noises (7,(t)n;(t')) are given by 


(ni(t) =| (Big — 0,5) (87,51 — gr NIEŻ ik (t-t) (1.23) 


U=1 


I have used the Langevin approach in the study of the stochastic elimination of fast 


variables [29] presented in Sec. II. 
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1.2.2 Master equation 
a) Reactive-only system 


Although a Langevin approach is a good start, it requires some assumptions on the size 
of the system and the nature of the fluctuations. The master equation offers a better 
description of a system at a mesoscopic scale [1, 14, 47]. We consider that the system has 


a given probability P({®},t) to be in a given state 
{8} = 4 Ny Nas. Nn} (1.24) 


at time t. The system state is described by the discrete numbers N; of molecules of each 
chemical species X;. During a finite time interval [t,t + 7], the probability of finding the 
system in a given state (f) evolves according to all possible reactions. The reactions 
are assumed to be Markov processes. Consequently, the probability P({8}, t+ T) for the 
system to be in state (f) at time t + 7 only depends on system state at time t and the 


transition rates or propensities between all states at time t and the state (f) at time t+7 


P({@},t+7) =P({8}, t) x Probability to remain in the state (f) 


+ | (19'),t) x Probability to jump from state (f') to state [91] (1.25) 
(9) 


The probability of leaving a given state is the sum of the probabilities to go from that state 
to any other one. The probability to remain in a given state is simply (1— Probability of 
leaving that state). If we note T({®} + (f'l) the transition rate from state (f) to state 
{P}, Eq. (1.25) can be written as 


P({®},t+7) — PUS}, t) = Y [P{#’}, t) x TUD) + 19)) — PUSY, t) x TUD) > 49'))] 
(2) 
(1.26) 


Intuitively, the probability of leaving the current state (f) is the probability that one 
reaction occurs in the time interval (t,t 4-7), as in Eq. (1.14). Therefore, we use the same 


assumption as in Eq. (1.15) but with discrete numbers of molecules 
S TUS} + {P} = Ler = 2 kj II Aa; T (1.27) 
(9) R, i 

where we made explicit the number of possible pairs of molecules 


N; Ni! 
Aas, = TF 1.28 
Oj (Ny = ay)! ( ) 
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Finally, assuming that the time 7 is sufficiently small to be considered as an infinitesimal 
time interval dt and using Eqs. (1.26) and (1.27), we write the master equation for a 


homogeneous reactive system 


m 


PUS}, t) > II Aa, AID pr), t) — II Adi PUB (129) 
il-1 


where the state (DP) is the state that evolves into the state (PD) after one elementary 


reaction j. 


b) With diffusion 


In order to take diffusion into account in a master equation, we introduce at least one space 
dimension and the transition rates describing the diffusive processes. For anv reaction- 


diffusion system, the master equation can be divided into two terms 


P({8}, t) = APP), leccion + OPD), t) |diffusion (1.30) 


where the first term accounts for the reactive processes described in Eq. (1.29) and the 
second term accounts for the diffusive processes. For the sake of simplicity, we assume that 
there are N molecules of only one chemical species X. The system is a one-dimensional 


L 
array of length L, divided into £ boxes of length Ax = 7 


o © © 
© © © e? O o 
A O o © © o © © 
ojo ojo o (e) © | 
0 Ar 2Azx L 


where each box is considered homogeneous. We use periodic boundary conditions. The 
number of molecules in each box i is denoted by N(i,t). The state of the system is given 
by 


(5) = {N(1,t), N(2,8),..., N(i,t), ..., N(£,t)) (1.31) 


Similarly to Eq. (1.26), the probability to leave the state (f) during a time interval |t, t+7] 
is equal to the probability that a single molecule of the system jumps from its current 
box to a neighboring one (1, 14]. Assuming that boxes are independent, the probability 


to leave the state (f) can be written as the sum of all probabilities in each box for one 


http://rcin.org.pl 


10 Chapter I. Methods 


molecule to jump 


Prob. to leave state {8} = X` T({8} > (9')) (1.32) 
fe! 
L 
— 5 Prob. that one molecule leaves box i (1.33) 
i=1 


The probability to remain in the state (f) is then simply (l-Probabilitv to leave state 
{®}). Hence, we only need to determine the probability that a molecule leaves a box i. 


This probability is given by the propensity d; that one molecule leaves the box i 


d;dt =(Probability that the molecule jumps to the left (1.34) 
+Probability that the molecule jumps to the right) (1.35) 
x Number of molecules in the box i (1.36) 


According to kinetic theory, the probability for a molecule to go left or right is proportional 
to the square root of temperature [14]. The Einstein relation gives the relation between 
the macroscopic diffusion coefficient of Eq. (1.3) and temperature. Typically, in large 


systems, the diffusive flux between two cells is related to the propensity d;. We have 


2D 
d; ~ 


i= Nit) (1.37) 


where D is the diffusion coefficient of the chemical species X. Using this approximation, 


we write the diffusion term of Eq. (1.30) for a single species in the form 


£ 
A PLP), t)|diffusion = 23 VIN; + 1) [PANG iz 1,t) sz 1, N(ż,t) + 1} F P({N(i, t) + 1, NG Sa 1,t) m 1)] 
i=1 


— 2N,P(45), t) (1.38) 


where only the number of molecules differing from N(i,t) in cell ¿ are made explicit. 

I have used a master equation approach in the study of the stochastic elimination of a 
fast variable [29] presented in Sec. II and in the analysis of fluctuation effects of an FKPP 
wave front [44] presented in Sec. IV.2.2. 


13 Concentrated solution 


Usually, in solutions, solvent particles are nonreactive and in large excess compared to 
reactants. Hereafter, a solution is said highly diluted when the local solvent concentra- 


tion is arguably constant. Considering the high-dilution limit ensures that the diffusion 
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coefficient of all reactants and the concentration of the solvent are constant. Hence, in 
the case of a reactive solvent, its concentration can be included into the rate constants of 
the corresponding reactions. However, small systems may be crowded, which invalidates 
the approximation of a large amount of solvent. The study of the effects of a deviation 


from the high-dilution limit in small systems is therefore necessary. 


1.3.1 High-dilution limit 


We present the effect of the deviation from the high-dilution limit on the reaction-diffusion 


dynamics using a simple example. We consider the Verhulst model [48] 
A+8—>28 (1.39) 


where A is a reactive species and S a reactive solvent. In the absence of hypothesis on 
the relative orders of magnitude of the concentrations cą and cg of species Å and S, the 


rate equation for A is 
dica = —kcącg (1.40) 


The total concentration c*** of chemical species 


ct = cą + eg (1.41) 
is constant according to Eq. (1.39). Therefore, we can eliminate Eq. (1.40)'s dependency 


on cg 


ctot 


dica = kc ca (i sA ) (1.42) 


and recognize the logistic function. The solution of Eq. (1.42) is 


ini — krtot 
dere kent 


ctot 4 ere SE 1) 


cą(t) = (1.43) 


where en is the initial value of ca. 

However, if the reaction scheme involves more steps and is more complex, the associated 
rate equations may not be solvable. In this case, it is possible to assume that the concen- 
tration of solvent cg is much larger than the concentration of reactants. Consequently, 


cg is considered constant and the reaction step becomes a lower-order reaction 


a 98 (1.44) 
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where the modified rate constant kg is given by 
kg = keg (1.45) 


This approximation is what we call the high-dilution limit. Our aim is to characterize the 


deviation from the high-dilution limit. Thus, we define the parameter 


bal c A 
i z S 
jz% (1.46) 


that represents the ratio of the sum of the concentrations c? of all non-solvent reactants 
evaluated at a steady state and the total concentration. When 6 tends to 0, the solvent 
is in large excess with respect to the other species and the high-dilution approximation 
is valid. On the contrary, when ô increases, the deviation from the high-dilution limit 
increases, and the previous assumption fails. This parameter is defined at a homoge- 
neous stationary state such that it remains constant even if the system exhibits spatial 
structures such as Turing patterns. If the system exhibits multiple stationary states, the 


parameter is defined according to the stable homogeneous stationary state of interest. 


1.3.2 Modified Fick’s law in a concentrated solution 


According to Fick’s law applied to a dilute system, the diffusive flux of one species is 
proportional to the gradient of its concentration. The associated diffusion coefficient is 
derived from kinetic theory according to the characteristics of the species. In a concen- 
trated solution, the concentration of one species has an impact on the diffusion of another 
species and cross-diffusion terms must be considered. We consider a concentrated solution 
of species A and B. The solvent S is still considered in excess but not in great excess, so 


that we are confronted with a ternary mixture of A, B, and S particles [36]. 


In a highly diluted solution, the center of mass of the solvent and the center of mass 
of the system are typically the same. However, we expect that, in a concentrated system, 
the two centers of mass are different. The idea is to exploit the frame of the solvent in 
which Fick’s law takes a simpler form [26, 36]. The flux of species X=A,B in the frame 


of the solvent is 
jx = cx(ux — us) (1.47) 


where uy is the velocity of the center of mass of species X in the frame of the system. 
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By definition, the flux of the solvent vanishes in the frame of the solvent 
j2=0 (Las) 


The fluxes in the frame of the solvent can be expressed in terms of the fluxes in the frame 


of the system using Eq. (1.47) 


is: = cx(ux —u+ u — us) (1.49) 
å Cx. 
= jx — js (1.50) 
cs 


which gives, using the law of conservation c'% = cą + cg + cs and Eq. (1.48) 


å CA |.S CA .S 

JA = (i 2) JA totJB (1.51) 
c c 

è CB „Ss CB V.S 

JB = tot JA 53 l= s) JB (1.52) 


Next step consists in using Fick's law relating the fluxes in the frame of the solvent and 
the concentration gradients. Within the framework of linear irreversible thermodynamics, 
the entropy production per unit mass due to isothermal diffusion is given by 

o=5 LY ix'(-Vrhx) (1.53) 


T X=A,B,S 


where T is the temperature, Vr the spatial gradient at constant temperature, and ux 
the chemical potential of species X. Assuming that the deviation from the high-dilution 
limit remains sufficiently small, the chemical potential for a given species is the same as 


in an ideal solution 
Cx 
ux = på + kgT log i (1.54) 


where på is the standard chemical potential of species X and kp is the Boltzmann constant. 


At constant pressure and temperature, the Gibbs-Duhem equation states that 


Y cx'(-Vrux)=0 (1.55) 
X=A,B,S 
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Using Eqs. (1.48) and (1.55), we write Eq. (1.53) in the form 
cx(ux — us + us — u) :(-Vrux) (1.56) 


jx-(-Vrux) (1.57) 


By introducing phenomenological coefficients Qxy to write linear relationships between 


fluxes and forces 


ix = Qxy(-Vrhv) (1.58) 
Y=A,B 


and using Eq. (1.54), we obtain the Fick's law in the frame of the solvent 


jx = 0) Divl-Ver) (1.59) 
Y=A,B 


where DŠ% are diffusion coefficients. However, even in a concentrated solution the no- 
tion of solvent keeps some relevance. Although the variation of the concentration of the 
solvent cannot be ignored, the concentrations of A and B are significantly smaller than 
the concentration of the solvent S. In these conditions, the vast majority of the binary 
collisions involve at least one S solvent particle. Hence, diffusion of species X=A,B is 
mainly imposed by the collisions between X and the solvent S while the impact of the 
collisions between A and B is negligible. We therefore admit that D$- is negligible for 
XY and we denote D3. y by Dż, for X=A,B. Consequently, Eq. (1.59) becomes 


j% = Dil-Vex) with X=A,B (1.60) 


so that the flux of X in the frame of the solvent only depends on the concentration of X. 
Finally, the modified Fick's law in the frame of the system that accounts for the deviation 


from the high-dilution limit in a ternary mixture is given by 


A CA CA 

JA == (i = 5) DSVca + aot DBVCB (1.61) 
A CB ns CB S 

JB = ¿tot DAVCA = (i = z) DBVcB (1.62) 


(1.63) 
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which can also be written under the matrix form 


PĘ (1-34) fa ta aq) (1.64) 


: c c 
JB =B |= SB -DÌ V CB 
(tot (tot 


I have used the modified Fick's law to study the effect of the deviation from the 
high-dilution limit in Sec. IV.2. 


1.4 Numerical methods 


1.4.1 Gillespie algorithm 


In the general case, the master equation presented in Eq. (1.29) is not solvable. Neverthe- 
less, it is possible to generate a stochastic trajectory using the Gillespie algorithm [49]. 
Gillespie uses a kinetic Monte Carlo procedure to directly simulate the reaction and dif- 
fusion processes and solve the master equation. The general formulation of the algorithm 


contains several steps, which are presented hereafter. 


a) Reactive-only system 


We consider the reaction scheme presented in Eq. (I.1) and the associated master equation 
given in Eq. (1.29). Each reaction step j has a certain probability to occur within a 
random time interval |, t + 7]. During this time interval, the propensity p; that the next 
elementary reaction is a reaction j is given by Eq. (1.14), with the same assumptions as 
for Eq. (1.27) 


n 
N; 
p; = k; Il 47, (1.65) 
il—1 


Thus, the propensitv po that anv reaction occurs is the sum of all propensities p; of the 


m steps 
m 
Po = Y pi (1.66) 
j=l 


If the time interval [t,t + 7] is split into infinitesimal intervals dt 


t EFT 


dt 


then the probability Po(T) that no reaction occurs in the time interval (t,t + 7] is equal 
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to 


= 
P>(r) (1 — podt) dt = e" Por (1.67) 


= lim 
dt—0 
Therefore, the probability Po that at least one reaction occurs is 
P) = 1 — e” (1.68) 
implying that the probability distribution of random reaction time 7 is 
p(t) = dr Py = poc” Po" (1.69) 
The random reaction time 7 is exponentially distributed with mean 
(== (1.70) 


The first step of the simulation is the initialization of all numbers of molecules, rate 
constants of the reaction, and the random number generators. Then, in the second step, 
we generate a random time interval according to Eq. (1.69) and select a random elementary 
reaction proportionally to its propensity using Eq. (1.66). In the third step, we update the 
number of molecules according to the reaction that occurred and increase the time step 
by the randomly generated reaction time. Finally, we go back to the second step where 
we generate a new random reaction time and a new elementary reaction. Eventually, the 
simulation stops when the number of reactants has reached zero or the simulation time 


has run out. 


b) With diffusion 


When diffusion processes are involved, the Gillespie algorithm can be easily adapted. We 
now consider the system described by Eq. (1.30). During a random time interval [t, t +7], 


the propensity d* that a molecule X leaves the box i is given by Eq. (1.37) 


2Dx 
2 


di = dé le + di right = Nx (i,t) (1.71) 


T 


where dZ ett and de |right denote the propensity that a molecule X in the box i jumps 
to the left or to the right, respectively. A priori, these two propensities are equal. The 


propensity do that any molecule leaves its box is 


l 
do = VE (1.72) 


X i=l 
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which can be added to the propensity given by Eq. (1.66). Then, we adapt the second 
step of the algorithm such that a random elementary process (reaction or diffusion) is 
selected proportionaly to its propensity. 


I have used Gillespie algorithm in the study of the elimination of a fast variable [ ] 
presented in Sec. II and the study of an FKPP wave front in a concentrated system [ | 
presented in Sec. IV.2.2. 


1.4.2 Direct Simulation Monte Carlo method 
a) Concept 


The Direct Simulation Monte Carlo method (DSMC), developed by Graeme Bird in 
the 60’s, is a numerical method used to compute molecular gas flows in aerodynam- 
ics{ , , |. It has been successfully extended to include reactive mechanisms and can 
be used to simulate highly diluted solutions. The method relies on a kinetic Monte Carlo 
algorithm which generates stochastic trajectories of particles and amounts to a direct 
simulation of the Boltzmann equations including fluctuations. Particles are hard spheres 
of mass m and radius r with continuous positions and velocities. Initial positions of the 
particles, compatible with the macroscopic initial conditions for the concentrations, are 
randomly chosen. The initial velocities are sampled according to a Maxwellian distribu- 
tion with kgT = 1. During a time step, positions are updated according to the velocities. 
The main feature of DSMC is to treat collisions statistically. The space is discretized into 
cells of length Az, where particles are susceptible to collide only with particles inside the 
same cell. According to the collision integral of the Boltzmann equation, the probability 
of collision of two particles is proportional to their relative velocity [ |. The "No Time 
Counter" (NTC) algorithm [ , | derives an integer upper bound to the maximum num- 
ber of collisions to be performed during the time step. An acceptance-rejection method 
is then used to test whether a collision between two randomly chosen particles in a box is 
accepted or not. The collisions are considered elastic from the mechanical point of view 
and the final relative velocity of the colliding pair is determined according to isotropic 
scattering. 

During a collision, a chemical reaction may happen. The reaction occurs with a proba- 
bility proportional to the corresponding rate constant determined by a steric factor pa 
and an activation energy Ea. According to kinetic theory, in a binary mixture of A and 
B species, the collision frequency is given by 


8kpT 
THAB 


ZAB = CACBOAB’ (1.73) 
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MB 
ma + MB 
of the reactants. Only the fraction of encounters that has a relative kinetic energy greater 


where cap is the cross-section of the collision and tap = is the reduced mass 


than the activation energy Fa of the reaction reacts. Therefore, the rate of the reaction 
is 
EĄ 
SZ A 1.74 
VAB ABPABEXP ( a) ( ) 


If we compare this result to Eq. (1.2), where the rate of a binary reaction is vag = kjcacp, 


the expression of the rate constant as a result of kinetic theory is 


8kgT 75) 
kan = 1.75 
AB = PABOAB\| WE exp ( RT (1.75) 


In the simulations that we performed, it is assumed that the activation energy of all 


reactions is equal to 0 and that no reaction is endothermic or exothermic. Therefore, 
the temperature is constant and homogeneous in all simulations. The cross-section of the 


collisions between two molecules Å and B is given by 
JAB = (ra + TB)” (1.76) 


for all collisions considered. The procedure used to obtain desired diffusion coefficients is 
also derived from kinetic theory. In a binary mixture, the diffusion coefficients of both A 


and B particles are equal and given by 


3 kgT 


Dy =D = DĄ Ez 
3 B 8(ca + cB) (ra + rp)? \ mua 


(1.77) 


In this expression, the diffusion coefficient depends on the local concentration ca + cp. 
However, diffusion coefficients are assumed to be constant in space and time in Eq. (1.4). 
This hypothesis requires that no particle is created ex nihilo or destroyed, and that ex- 
changes with the exterior (such as reservoirs) do not radically change the local concen- 
tration. Therefore, we make sure that the concentration cą + cp is arguably constant in 
the simulation. Another condition for collisions to be correctly simulated is that the cell 


length Ax is smaller than the mean free path of colliding particles 


1 


(= I.78 
v2ctotr(r4 + rg)? ( ) 


In highly diluted solutions, the excess of solvent makes this condition hard to fulfill while 
keeping reasonable computation times. The condition can be relaxed if the mean gradients 


of concentration between two neighboring cells are typically smaller that the amplitude 
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of the concentration fluctuations. 
I have used the DSMC method to simulate species with different diffusion coefficients and 


Turing structures in a concentrated system [37]. The results are presented in Sec. II.3. 


b) Third-order reaction 


In the article [37] presented in Chapter III, we consider a ternary mixture of A, B, and S 
species where one step of the reactive mechanism is a third-order reaction. The difficulty 
was to adapt the DSMC algorithm to correctly simulate a third-order reaction. The issue 
has been solved as follows by M. Mareschal eż al. [51] and already used by the group [52]. 


The reaction 
A REST SĄ (1.79) 


is split into two second-order reactions 


ABL, AB (1.80) 
ABA PZ SĄ (1.81) 


where the first reaction occurs with the same rate constant k as the third-order reaction 
and the second reaction is supposed instantaneous. The simulations reproduce the third- 
order reaction given in Eq. (1.79) as follows. When a binary collision between a particle 
A and a particle B is accepted in a given spatial box, a complex AB is supposed to be 
formed. A third particle is randomly chosen in the same spatial box. If it is an A particle, 
the AB complex is immediately transformed into two A particles. Hence, according to 


Eq. (1.75) and the probability = of picking a particle A, the rate constant k obeys 
c 


4 : 1 
k= pa ale y nkgpl x — (1.82) 


(tot (tot 


with wap = Ea = pap = 1. Due to the proportionality to 1/c'%, the rate constant k 
depends on the deviation 6 from the high-dilution limit from Eq. (1.46). 


c) Diffusion in a ternary mixture 


Observing a Turing structure requires the simulation of sufficiently different diffusion coef- 
ficients for the activator and the inhibitor. As already mentioned, the diffusion coefficients 
of the two components are identical in a binary mixture. The solvent, introduced to study 
concentrated solutions, advantageously plays the role of the third kind of particles with 
respect to diffusion. The problem of tuning the diffusion coefficients in the simulation of 


a ternary mixture has already been addressed by the group [53] and I recall the main lines 
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of the method. Qualitatively, faster diffusion of the inhibitor is obtained by considering 
particles of the same mass but different diameters. In a ternary mixture of A, B, and 


S [54], the fluxes of matter ją, jg, and js are expressed as 


for X,Y, Z = A, B, R, where the three-component diffusion coefficients Dg, Dys, and 
Da g are derived from the two-component diffusion coefficients Dag, Das, and Dgs from 
Eq. (1.77) 


cz(Dxz — Dxy) 


Dyy = Dxy Il- 
oo cx Dyz +cyDxz +czDxy 


(1.84) 


The presence of the concentrations in the expression of the diffusion coefficients may imply 


non desired space-dependent diffusivities. However, considering 


cgDyp > caDps + cBD' 49 (1.85) 
Cs > ca (1.86) 
Cs > CB (1.87) 
and using Eq. (1.84) leads to 
Dan = Das = Das (1.88) 
Dra ~ Dips ~ Dgs (1.89) 


which, combined with Eq. (1.83) for ¢°t = cą + cg + cg constant, gives 


ja = —Das0;CcA (1.90) 
jp = —DBsózcB (1.91) 


Thus, when the solvent is sufficiently in excess, the diffusion coefficients of Å and B 
species are the same as in a binary mixture of Å and S particles and B and S particles, 
respectively. Intuitively, in a system with S particles in excess, the proportion of collisions 
involving no S particles is negligible. Hence, we assume that the diffusive mechanism of 
A and B particles is dominated by their interaction with the S particles. For the sake of 


simplicity, we write 


Da = Das (1.92) 
Dp = Ds (1.93) 
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In Eq. (1.77), the diffusion coefficient depends on the local concentration of the two 
species involved in the diffusive process. In Chapter. III, we study Turing structures with 
spatial oscillations of the concentrations, which undermines the hypothesis of constant 


diffusion coefficients. In order to derive appropriate diameters for the different types of 


D 
particles [53], we consider spatial averages of the concentrations. Then, writing d = A 
A 
and ma = mg, and using Eq. (1.77), we get 
l—vd 
PE + (1- vd)rs (1.94) 
vd’ 
d 
where d = d(es + ca) ~d. As rg must be positive, we obtain the condition 
CR + CB 
TĄ > (Vd! —1)rg (1.95) 
If this last condition is satisfied with a small margin, then 
TB KTA (1.96) 


and Dp > Da as desired. 
Finally, introducing Eq. (1.95) into Eq. (1.85), we obtain the condition on the concentration 


of solvent S for the proper simulation of the diffusion coefficients 


cs > V(dea +cp)(ca + cp) (i — 75) (1.97) 
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Chapter II 


Stochastic approach to the steady-state 


approximation 


The usual adiabatic elimination often encountered in chemistry is the steady-state ap- 
proximation, consisting in eliminating a fast concentration. In a complex mechanism, 
identifying a fast concentration is not straightforward [55, 56, 57, 58, 59, 60]. A linear 
analysis may be locally performed and requires linearizing the rate equations, computing 
the eigenvalues, and using the change-of-basis matrix to relate the concentrations and the 
eigenmodes. The relationships between the eigenvalues and the rate constants may not 
be trivial and the knowledge of the rate constants is not always sufficient to identify a 


fast variable at first glance. 


11.1 Context 


We give an example of steady-state approximation in a simple case involving two elemen- 
tary steps with rate constants of different orders of magnitude, sufficient to generate a 
fast concentration. The example also illustrates how third-order steps may be obtained 
by reduction of a mechanism containing second-order steps. 

We consider a reaction scheme involving two elementary steps and four species A, B, 
C, and D 


A+B 25 C 
A+C 0 (ILL) 
The first step is supposed to be much slower than the second one, i.e. the parameter € 


obeys £ < 1. The quantities cą + cc + €p and c4 — 2cp — cc are conserved and dynamics 


involves two independent variables. It is however simpler to keep the three variables ca, 
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cp, and cc. The rate equations are given by 


dca = —ekcacp — kcąco (11.2) 
d;cp = ekcacp (11.3) 
dcc = ekcacp — kcacc (11.4) 


Introducing a new time scale 7 = et, expanding the concentrations into power series of e, 


and using Eq. (11.4) at zeroth order, we deduce that W = 0. Consequently, we have 


del) = ke (co) — c) =4) (11.5) 


which implies c = cO, Finally, Eqs. (11.2) and (11.3) lead to 


d,ch) = —2kc4(0)ep(0) (11.6) 
d,ch) = —kca(0)cg(0) (11.7) 


which corresponds to the rate equations associated with a third-order reaction 
LAB 59D (11.8) 
The reaction between nitrogen monoxide and chlorine 
2NO + Ch — 2 NOCI (11.9) 


illustrates the reduced mechanism given in Eq. (11.8). The mechanism proposed in Eq. 
(11.1) does not rely on any chemical considerations and is only one possible two-step 


mechanism compatible with the reaction between nitrogen monoxide and chlorine. 


11.2 Summary of the results 


In the previous example, the dynamics of the system is described at a macroscopic scale. 
The nonlinearities of the deterministic dynamics and the large fluctuations reached in 
small systems are known to interact [61, 62, 63] and make the elimination of fast concen- 
trations a nontrivial problem in stochastic systems [58, 64, 65]. In biological experiments, 
Fluorescence Correlation Spectroscopy (FCS) [66] is widely used to study the dynamics 
of labeled species, for example, to evaluate rate constants when the reaction scheme is 
supposed to be known. The interpretation of the results requires the comparison of the 
experimental data with analytical expressions of the correlations of concentration fluc- 


tuations for the reaction scheme of interest. However, reaction schemes in biology often 
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involve hundreds of steps and the computation of the correlations is not tractable before 
a reduction of the mechanism. Defining the scope of fast concentration elimination in 
stochastic nonlinear dynamics is therefore essential to interpret FCS results. Similarly, 
the reduction of a mechanism may lead to nonpolynomial nonlinearities, as exemplified 
by the reduced Michaelis-Menten scheme. The conclusions that could be derived from 
the stochastic analysis of the reduced Michaelis-Menten model may differ from the direct 
analysis of the complete scheme [67]. In order to analyze the consequences of mecha- 
nism reduction, I consider a minimal chemical model involving two species of variable 
concentrations capable of evolving into a Turing pattern. Then, the two-variable model 
is assumed to be the reduction of two different three-variable models. The problem is 
to determine if the correlations of fluctuations in the three-variable models are correctly 
predicted by the two-variable model, in the limit where the reduction of determinis- 
tic dynamics is valid. I developed an analytical approach based on chemical Langevin 
equations linearized around the steady state of interest as presented in Sec. 1.2.1. Fol- 
lowing the method applied in references [61, 62, 63] to characterize the asymmetry of 
time cross-correlations in far-from-equilibrium systems, I determined the expressions of 
the correlations of concentration fluctuations in the two- and three-variable models. In 
parallel, I performed simulations of the corresponding master equations (see Sec. 1.2.2) 
using Gillespie algorithm according to the procedure given in Sec. 1.4.1. The weaknesses of 
the Langevin approach in the description of the internal fluctuations in a nonlinear chem- 
ical system have been highlighted [68]. The master equation approach has proven that 
the two-variable model does not correctly predict the fluctuations in the three-variable 
systems. We concluded that the coupling between the fluctuations and the nonlineari- 
ties of deterministic dynamics makes the use of the steady-state approximation delicate 
when the studied system requires a good control of the fluctuations. The predictions 
of the correlations based on a reduced mechanism must be considered with special care 
when preventing hazards in explosive phenomena, modeling pattern formation in biol- 
ogy, or dealing with small systems in which variances of fluctuations are detected as in 


fluorescence correlation spectroscopy [65, 69]. 


11.3 Publication 


The results are published in the article “Elimination of fast variables in stochastic non- 
linear kinetics”, G. Morgado, B. Nowakowski, and A. Lemarchand, Phys. Chem. Chem. 
Phys., 22, 20801-20814 (2020) [29]. 
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Chapter III 


Turing patterns 


In 1952, Alan Turing presents one of his later most famous works [5]. He proposes 
a mathematical analysis of reaction-diffusion systems that may exhibit periodic spatial 
structures from an initially homogeneous state. The discovery of time oscillations in 
a chemical reaction was made in the 50s by the Russian biochemist B. Belousov, while 
looking for a non-organic analog to the Krebs cycle [70]. A young chemist, A. Zhabotinsky, 
became famous in the 60s for the work he devoted to this reaction with the observation 
of spatial structures and chemical wave fronts [71]. Nowadays, the tradition associates 
their two names to designate this complex reaction. Experimental evidence of Turing 
type nonequilibrium chemical patterns has been provided in 1990 by the group of De 
Kepper [72]. 


111.1 Emergence of a Turing pattern 


Turing patterns (or Turing structures) only appear in out-of-equilibrium systems, there- 
fore relying on the outside environment to maintain the pattern. Initial conditions and 
parameters such as the kinetic rate constants and the diffusion coefficients of the reac- 
tive species are also crucial for the existence and the shape of the pattern. The patterns 
appear when an inhomogeneous perturbation arises in a homogeneous stationary state. 
One of the simplest systems that exhibits Turing patterns is the Gray-Scott model [73], 


which involves two reactive species A and B in a third-order autocatalytic reaction 


EF (ILL) 
2A+B 34 (111.2) 
B == Ro (IIL3) 
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where Ra and Rp are reservoirs of particles A and B, respectively. The associated 


reaction-diffusion equations are 


Oca = —kica + kącącB + DąOŻCA (111.4) 
cB = —kocĄcB — kscp + k-3 + Dpóźcp (IIL5) 


where k_3 = k*3cr, with cr, constant. The system possesses three stationary states. 
The steady state 


(111.6) 
pe (III.7) 
corresponds to the absence of A particles and is obviously stable according to the chemical 


scheme given in Eqs. (111.1-111.3). The two other states are derived from Eqs. (III.4) 
and (111.5) 


cj = (IIL.8) 


LI 
TT 
II 
o 
w 
l 
GON 
fe 
a 
mH 
jiġ 
W 
+H 
> 
IN 
w 
l 
Aa 
> 
ar ED 
a 
w 


= (111.9) 


By performing a linear stability analysis, as presented in Sec. 1.1.3, we show that the 
state (cą, c$) is stable towards homogeneous perturbations whereas (cą, Cp) is unstable. 
Knowing that Turing patterns emerge from inhomogeneous perturbations, we consider 
the evolution of a local inhomogeneous perturbation dc = (dca,dcp) around the state 
(ci, ch). According to Eqs. (III.4) and (III.5) and in the framework of a linear approach, 


the Fourier transform of the perturbation 
je I dzóc - e 1 (III.10) 
VIT 
obevs 


DIĠCA = —I40G4 + ka (2chcHóCA + ch den) —g'Dadca (111.11) 


OER = —ka (264c56CA + CĄÓG) — ksden — q*DnóGR (111.12) 
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CB, Stable 
(ch, 2) 
ko Jnstaple 
TA TA 
w" K 
` 
kz = Stable 
A. (ch, cp) 
Ca 
5 © 
~ 
v. 
v. 
0 k-s fa CA 
2ki ki 
Fig. III.1 Stationary states of the model given in Eqs. (III.1-III.3) 
in the space of concentrations (ca, cp) 
from which we obtain the stabilitv matrix 
- 0D 
M = lis dec a (111.13) 
mig mas — q DB 
where 
miu = —k1 + 2k2ch ch (IIL.14) 
m21 = hack (111.15) 
mig = —Żkącąch (111.16) 
mas = —kąch” = k3 (111.17) 
The two eigenvalues of the stability matrix are 
_ 1 2 A. 2 2 
Mz = 5 [Mi + mo — 9 (Da + Dg) = y (m — m22) — (Da — DB)] + 4m21m32 
(111.18) 
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The steady state (di ch) is unstable if the real part of one eigenvalue is positive. We 


look for the Fourier mode qmax that maximizes the largest eigenvalue uy 


= Hii De [= 
TA ma — M22 AT DB m12M21 (111.19) 
Da-Dg  Da-DgYV DAĄDB 


If the value p14 (max) is positive, then the inhomogeneous perturbation grows and a Turing 


pattern emerges. The pattern consists of sinusoidal oscillations of Å and B concentrations 


with the wavelength 


27 


A= (111.20) 


dmax 


Equation (III.19) is valid for any reaction mechanism involving two chemical species sus- 
ceptible to evolve into a Turing pattern. We can therefore express general conditions for 
the emergence of Turing patterns. First, we remark that if the two diffusion coefficients 
are equal, i.e. Da = Dp, then Eq. (111.19) diverges and no Turing pattern can emerge. 
Second, the wavelength depends only on the rate constants and the diffusion coefficients, 
and not on the boundary conditions of the system. This makes the Turing pattern inde- 


pendent of the size of the system, at least at the macroscopic scale [74, 75, 76]. 


111.2 Termination mechanism of Turing patterns 


11.2.1 Context 


Models of periodic spatial patterns usually involve infinite systems or periodic boundary 
conditions (77, 78, 79]. However, in morphogenesis, the question of the termination of 
a structure arises [80, 81]. Specifically, the spine of the vertebrates ends with smaller 
vertebrae. In the framework of Turing patterns, this phenomenon implies both a decrease 
of the amplitude of the spatial oscillations and a decrease of the wavelength. Deciphering 
the mechanisms controlling the termination of the spine in an embryo is far beyond the 
scope of this work. Our aim is to propose a possible mechanism at the macroscopic scale, 
inspired by biological structures and compatible within a chemical engineering context. 
The boundary conditions chosen by the group in 2016 [36] to model the growth of the 
spine are well adapted to the design of an artificial spatial structure. In the model of 
reference [36], the Turing pattern develops behind a wave front propagating from left to 
right. Neumann or second-type boundary conditions are chosen at the left boundary at 
which the derivative of the concentration with respect to the spatial coordinate x vanishes. 
The emerging Turing pattern will therefore possess an extremum at the left boundary. A 


free boundary is imposed at the right, i.e. the system grows freely in this direction, so 
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that the propagation of the wave front to the right is not perturbed. The right part of 
the Turing pattern, located behind the front, is free and not affected by the boundaries. 


11.2.2 Summary of the results 


As explained in Sec. III.1, the stability of the Turing pattern is related to the eigenvalues 
of the stability matrix around the homogeneous steady state. A real, positive eigenvalue 
y corresponds to an unstable homogeneous steady state and a stable Turing structure. 
As uy tends to 0, the amplitude of the spatial oscillations decreases. Destabilization is 
reached when w+ = 0. Using Eqs. (111.18) and (111.20), I performed a systematic analysis 
of the effect of all rate constants and diffusion coefficients of the model given in Eqs. (III.4) 
and (III.5) on the eigenvalue u; and the wavelength of the structure A. Interestingly, a 
monotonous variation of almost any of the dynamical parameters leads to the desired 
behavior. In particular, either the increase or the decrease of a rate constant leads to the 
simultaneous loss of stability of the structure and the decrease of its wavelength. Only 
the variation of the diffusion coefficient Da of the activator causes anti-correlated results, 
i.e. a decrease of the oscillation amplitude and an increase of the wavelength. For a given 
chemical system, locally varying a rate constant or the diffusion coefficient Dg of the 
inhibitor is not straightforward in the framework of chemical engineering. For an easy 
implementation, I suggest to impose an appropriate spatial profile for the concentration 
of the reservoir Rp, resulting in the increase of the effective rate constant k_3Rg of the 


process given in Eq. (111.3) and the desired termination of the structure. 


111.2.3 Publication 


The results are published in the article “Termination mechanism of Turing patterns in 
growing systems”, G. Morgado, L. Signon, B. Nowakowski, and A. Lemarchand, Acta 
Phys. Pol. B, 50, 1369 (2019) [39]. 
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The question of the termination of a periodic spatial structure of Turing- 
type in a growing system is addressed in a chemical engineering perspective 
and a biomimetic approach. The effects of the dynamical parameters on the 
stability and the wavelength of the structure are analytically studied and 
used to propose experimental conditions for which a Turing pattern stops 
by itself with a decreasing wavelength. The proposed mechanism is suc- 
cessfully checked by the numerical integration of the equations governing 
the dynamics of the activator and the inhibitor. We conclude that a local 
increase of the concentration of the reservoir which controls the injection 
rate of the inhibitor into the system can be used to achieve the appropriate 
termination of a Turing pattern. 


DOI:10.5506 /APhysPolB.50.1369 


1. Introduction 


During embryonic development, segmented structures of the body such 
as the spine and the digits are formed by the production of repeated ele- 
ments. Since the seminal work of Turing |1] accounting for the formation 
of biological pattern in the framework of reaction-diffusion models, exper- 
imental evidences of Turing structures have been given in chemistry [2-4] 
and biology (5, 6]. Recent years have witnessed a growing number of pa- 
pers where reaction—diffusion principles are proposed to model the forma- 
tion of biological periodic spatial structures [7-13]. Following Turing, a 
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two-component chemical system composed of an autocatalytically-produced 
activator by consumption of an inhibitor that diffuses faster may produce 
periodic spatial structures such as stripes in one-dimensional (1D) systems 
and hexagons in 2D. In other words, a Turing pattern emerges by local 
self-activation and lateral inhibition [14]. The concepts developed to model 
living systems provide a source of inspiration in chemical engineering (15-22). 
However, standard models of Turing patterns generate structures in infinite 
systems and the question of the termination of a striped structure in a finite 
system arises in a perspective of biomimetism in material science. Specifi- 
cally, it is desirable to find experimentally achievable conditions creating a 
finite-size structure, whose growth stops by itself with decreasing oscillation 
amplitude and respects the decrease of the wavelength during the termina- 
tion process. To this aim, we use an as simple as possible reaction—diffusion 
model [23| admitting a Turing structure developing behind a propagating 
wave front and examine the effect of all parameters on both the stability 
and the wavelength of the structure [5, 22]. We already used a stochastic 
approach to a Turing pattern [23] and showed that, contrary to intuition, the 
internal fluctuations may have a constructive effect and stabilize the struc- 
ture outside the domain of stability predicted by a deterministic description. 
Here, we adopt a macroscopic approach. Our goal is to select suitable condi- 
tions from this systematic approach and to propose termination mechanisms 
compatible with processing in chemical engineering. 

The paper is organized as follows. In Section 2, a reaction-diffusion 
model involving local activation and long-range inhibition is presented. An 
analytical stability condition and the wavelength expression of the Turing 
pattern are given. The influence of the parameters of the model on the 
stability and the wavelength of the pattern is studied in Section 3. The 
analysis of the results leads to the selection of a user-friendly termination 
mechanism in the framework of chemical engineering. The analytical predic- 
tions regarding stability and wavelength are compared to numerical results 
for the chosen mechanism. Section 4 contains conclusions. The possibil- 
ity that the different mechanisms exhibited could be found as termination 
scenarios in morphogenesis is raised. 
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2. Model 


We consider the following reaction mechanism inspired from the Schna- 
kenberg model [24] and the Gray-Scott model [25]: 


A EE Ru, (1) 


AB 2% 34, (2) 


B = Ba (3) 
k_3 

where R; and Ro are reservoirs. The concentrations Ri and Ro of the reser- 
voirs are assumed to be constant in time. The reaction given in Eq. (2) 
autocatalytically produces species A and consumes species B. Due to the 
ability of accelerating its own production, species A is called an activator 
whereas species B, which is consumed by the same process, is called an in- 
hibitor. The macroscopic dynamics of the system is governed by two partial 
differential equations [9, 23] 


OA A 

EA Da 

At ka + ko + A Ox? ) (4) 
OB a 07 B 

a k_3Ro — k3B — k2 A” B + DB (5) 


for the concentrations A and B of the activator and the inhibitor supposed 
to have different diffusion coefficients Da and Dg. For appropriate rate 
constant values, such that 


A = (k_3R2)? — 4kikg/ko > 0, (6) 
the system admits two steady states (Ag = 0, Bo = k-3R2/k3) and 
k_3R2+ VA 
Ar = == 7 
T 2k1 , ( ) 
k-3Rə — V A 
By = ER? VA (8) 
2k3 


that are stable with respect to homogeneous perturbations. The index T 
stands for Turing. A linear stability analysis of Eqs. (4) and (5) reveals that 
the steady state (Ar, Br) can be destabilized by inhomogeneous perturba- 
tions (3, 5, 9, 23]. The Fourier transforms Aq(t) = [E A(z,tje “dz and 
B,(t) = [E B(u,tje "dx, where q is the Fourier mode, are introduced. 
In the Fourier space, the linear stability operator M is given by 


_ (ki — Dad? kin At, 
M = ( TA EA Dee (9) 
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The eigenvalues of the matrix M obey the characteristic equation u? + au + 
B= 0, with a = ky — 2528 Ar- (Da +Dp)q? and B = 2k? Ar/Br — (kı — 
Daq’)(k_3R2/By + Dgq?). The Turing structure develops if the largest 
eigenvalue 


1 kok 
(r AR: gr (Da + Da) 


K=J ki 


kok-3Ro A : 
d ki 4 4, AT T (DB — Da) a Sky ko Af (10) 
1 


is real and positive [3, 5]. Indeed, a system of differential equations, lin- 
earized around a homogeneous steady state, is easily solved by diagonalizing 
the linear operator. Then, the solution is a linear combination of eigenmodes 
which exponentially depend on time according to the corresponding eigen- 
values. Å term associated with a real, positive eigenvalue grows in time, 
leading to trajectories in the concentration space that move away from the 
fixed point [5]. In the studied system, the destabilization of the steady state 
occurs in favor of a Turing pattern. Equation (10) imposes conditions on 
the parameter values. In particular, the diffusion coefficient Dp of the in- 
hibitor B must be sufficiently larger than the diffusion coefficient Da of the 
activator A: The destabilization of the homogeneous steady state (Ar, Br) 
requires local self-activation, ensured by the autocatalytic production of the 
activator through the reaction given in Eq. (2), as well as long-range inhibi- 
tion, due to the larger diffusion coefficient of the inhibitor. The mode qmax; 
which maximizes the eigenvalue u, is the most unstable Fourier mode 


= > + Dp)v/2kik2DA/Dp— ki — kok-3RoAr/ki (11) 


Dp — DĄ 


In order to account for the termination of the Turing pattern in a growing 
system, including the fact that the structure ends with a gradually shorter 
spatial oscillation, we need to find conditions for which the structure tends 
to lose its stability while its wavelength decreases. The wavelength of the 
periodic structure is given by 


27 


= (12) 


dm ax 


The Turing structure becomes unstable as the value of the largest eigenvalue 
vanishes for the mode qmax associated with the maximum of u 
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1 kok aR 
1 (r RAE es JE 


Umax = 2 ki 


kąk_3Ro zo ; 
+ kj + 4, ^r + (DB — DAGE ax. — 8k ko AF (13) 
1 


with qmax given in Eq. (11). Figure 1 illustrates the behavior of Umax for 
parameter values associated with a stable Turing pattern with Umax > 0. It 
is also shown that it is sufficient to increase the value of k_3Ro to shift the 
curve (q?) down and lose the stability of the Turing pattern. 


1 T T 


0.5 F = 


0.5 tpai l | l [Ny 


Fig. 1. The largest eigenvalue u of the linear operator M versus square of Fourier 
mode q?. Solid line: k_3Ra = 8.76. Dashed line: k_3Ra = 10. Other parameter 
values: ki = 2.92, kə = 1, k = 2.19, Da = 1, Dg = 10. 


In the next section, we investigate the behavior of A and Umax as each 
parameter controlling dynamics varies. Specifically, we aim at identifying 
diffusion coefficients or rate constants whose variation leads both to a de- 
crease of the wavelength and a destabilization of the Turing structure, i.e. 
negative values for the maximum of the eigenvalue. 


3. Results 


The concentration Rə of the inhibitor reservoir is first assumed to be 
homogeneous. Figures 2 and 3 show the variations of the wavelength A 
and the maximum value max of the eigenvalue with one of the diffusion 
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Fig.2. Top: Scaled wavelength A/Az of the Turing pattern versus diffusion coeffi- 
cient Da of species A. Bottom: Maximum value max of the largest eigenvalue of 
the linear operator M versus Da. Parameter values: ki = 2.92, ko = 1, k3 = 2.19, 
k_3 Rə = 8.76, Dp = 10, Ax ="0.L; 
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Fig.3. Top: Scaled wavelength A/Az of the Turing pattern versus diffusion coeffi- 
cient Dg of species B. Bottom: Maximum value max of the largest eigenvalue of 
the linear operator M versus Dg. Parameter values: kj = 2.92, ko = 1, k3 = 2.19, 
k_3 Rə = 8.76, Da = ibs Ax = Dil, 
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coefficients, the other parameters being constant. The results are deduced 
from Eqs. (11) and (12) for A and Eq. (13) for max, the expressions of the 
steady state (År, Br) being given in Eqs. (7) and (8). To facilitate the 
comparison with the numerical integration of Eqs. (4) and (5) that will be 
performed in the following, the wavelength is given in a number of spatial 
cells of length Ax = 0.1. As shown in Fig. 2, the decrease in the maximum 
value Umax Of the eigenvalue as the diffusion coefficient Da of species A 
increases is accompanied by an increase of the wavelength A: The loss of 
stability of the Turing structure occurs with an increase of the spatial period. 
We conclude that a variation of the diffusion coefficient D 4 cannot be argued 
as a justification of the termination process. The behavior with respect to 
the diffusion coefficient Dg of species B is different. The simultaneous loss 
of stability of the structure and the decrease of the wavelength are observed 
in Fig. 3 as Dp decreases: The diffusion coefficient Dg of species B can be 
considered as a suitable parameter in the search for a termination model. 


Figures 4—7 show the variations of the wavelength A and the maximum 
value max Of the eigenvalue with rate constants. The variations of Å are 
given in a bounded interval of rate constant values, in which the Turing pat- 
tern is stable. At one of the endpoints of the interval, the eigenvalue [max 
vanishes and at the other endpoint, the condition of existence of the steady 
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Fig. 4. Top: Scaled wavelength A/Az of the Turing pattern versus rate constant kı 
of the chemical reaction given in Eq. (1). Bottom: Maximum value max of the 
largest eigenvalue of the linear operator M versus kj. Parameter values: kg = 1, 
k3 = 2.19, k_3Ra = 8.76, Da a 1, DB | 10, Ay €10, 
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Fig. 5. Top: Scaled wavelength A/Az of the Turing pattern versus rate constant ko 
of the chemical reaction given in Eq. (2). Bottom: Maximum value Lmax of the 
largest eigenvalue of the linear operator M versus ka. Parameter values: ki = 2.92, 
k3 = 2.19, k_3Ra = 8.76, Da = E Dg = 10, Ar = 0.1. 
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Fig.6. Top: Scaled wavelength A/Ax of the Turing pattern versus rate constant kg 
of the forward chemical reaction given in Eq. (3). Bottom: Maximum value [max 
of the largest eigenvalue of the linear operator M versus k3. Parameter values: 
ky = 2.92, ko = Il. k-3 Ro = 8.76, Da = 1, Dp = 10, Az =M 1 
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Fig.7. Top: Scaled wavelength A/Ax of the Turing pattern versus rate con- 
stant k_3R2 of the backward chemical reaction given in Eq. (3). Bottom: Maxi- 
mum value max Of the largest eigenvalue of the linear operator M versus k_3Ro. 
Parameter values: kj = 2.92, ka = 1, kg = 2.19, DA = 1, Dp = 10, Az = 0.1. 


state (Ar, Br) given in Eq. (6) is no longer satisfied. The two desired behav- 
iors, i.e. the decrease of both A and Umax, are observed as ki decreases, ko 
increases, kg decreases, and k_3Ro increases. For an assumed homogeneous 
concentration Ro of the reservoir, the variations of A and Umax with Ro are 
analogous to the variations with k_3R3. According to the chemical reaction 
given in Eq. (1), decreasing the rate constant ki tends to increase the con- 
centration of species A. Following Eq. (2), increasing the rate constant ka of 
the autocatalytic step tends to increase the concentration of species A and 
decrease the concentration of species B. This last result seems to be incon- 
sistent with the consequences drawn from the decrease in k3 or the increase 
in k_3Ro, which result in increasing the concentration of species B according 
to Eq. (3). However, we already stated that increasing B through soliciting 
the reservoir Ro results in consuming species B faster by the autocatalytic 
step given in Eq. (2) [9, 26]. In particular, we observed that introducing a 
local source of species B leads to the nonintuitive local decrease of B con- 
centration. Hence, all the variations of the rate constants that lead to a loss 
of stability of the Turing pattern are eventually associated with an increase 
of A concentration and a decrease of B concentration. 
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The diffusion coefficients and the rate constants characterize dynamics 
and are intrinsic to the reaction—diffusion system. Nevertheless, it is always 
possible to imagine spatial variations of the dynamical parameters. Well- 
chosen variations of the diffusion coefficient Dp of the inhibitor and each of 
the four rate constants of the chemical mechanism could be a priori used 
to build a termination model. In the framework of the application to devel- 
opmental biology, steric hindrance and molecular crowding in the tail of an 
embryo may be invoked to justify the decrease of the diffusion coefficients. 
In chemical engineering, a local increase of temperature could be used to 
induce a local increase of the rate constants. However, a local increase of 
molecular crowding or temperature is susceptible to simultaneously affect 
several dynamical parameters (13, 22, 27-32]. Whereas a decrease of Dg is 
desired to destabilize the Turing pattern, while decreasing its wavelength, a 
simultaneous decrease of DĄ would be detrimental. Similarly, an increase 
of ko and k_3 due to temperature increase could be satisfying but the joint 
decrease of kj and k3 could blur the effect on the Turing structure. The 
simplest way to imagine the control of a targeted parameter leading to the 
desired behavior is to impose well-chosen spatial variations of the reservoir 
concentration Ro. Indeed, the product k_3Ro plays the role of an apparent 
rate constant for the backward reaction given in Eq. (3) that can be fixed 
at will in chemical engineering in the case of a single dynamical system with 
uniquely defined intrinsic parameters. 

According to Fig. 7, increasing Ro tends to destabilize the Turing pattern 
and decrease its wavelength. We examine if the results deduced from a 
stability analysis can be used in a dynamical approach. The results of the 
numerical integration of Eqs. (4) and (5) for a homogeneous concentration 
Ro and a piecewise linear profile are given in Fig. 8. The initial condition is 
a step function between the steady state (Ar, Br) in the first 10 cells on the 
left and the steady state (Ag, Bo) in the remaining cells. The initial total 
number of cells is set at no = 610. Introducing the cell label i = x/Az, 
where Az is the cell length, and the discrete time s = t/At, where At is the 
time step, we choose 


A(i,s=0) = Ar, B(i,s = 0) = Br, for 1<i<10, (14) 
A(i,s =0) = Ao, B(s,s = 0) = Bo, for 11 <i< no. (15) 


To account for the growth of the system and simultaneously avoid boundary 
effects that may alter the wavelength of the structure [16], we impose a fixed 
boundary on the left and a free growing end on the right (9, 23, 26]. For 
parameter values for which the steady state (Ar, Br) is unstable with respect 
to inhomogeneous perturbations, a Turing pattern develops after the passage 
of a chemical wave front. More precisely, according to Eqs. (4) and (5) and 
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due to the no-flux boundary conditions applied on the left boundary, the 
concentrations in the first cell obey 


A(1,s+1) = A(1,s) — kjAtA(1, s) + ko At A(1, s)?B(1, s) 


+ Dara (As) = AN, (16) 
B(1,s+1) = B(1,s) + k_3RaAt — ksAtB(1, s) — ko At A(1, 5)? B(1, 8) 
Dr BL.) — B(1,5)) (17) 


so that both A(1, s) and B(1,s) are extremum of the Turing pattern in the 
first spatial cell į = 1. 
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Fig. 8. Spatial profiles deduced from the numerical integration of Eqs. (4) and (5) 
for kj = 2.92, ka = 1, kz = 2.19, k_3 = 8.76, Da = 1, Dp = 10, At = 10-4, 
tena/At = 2000000, Az = 0.1. Black dotted line: Imposed concentration Ro of 
the reservoir. (a) Homogeneous concentration Rə = 1, (b) Piecewise linear Ro 
profile. Gray/blue dashed line: Concentration of species A versus cell label x/Ax. 
Black/red solid line: Concentration of species B versus cell label z/Az. 


Spatial cells are added to the right end of the system at the front speed to 
counterbalance the progression of the wave front and mimic system growth: 
At all the discrete times s for which the concentration B(n — 600, s) of 
species B in the n— 600 cell becomes smaller than 0.99 Bo, the total number n 
of cells is increased by 1. Provided that the front propagates at a speed 
smaller than Ax/At, this protocol ensures that a layer of about 600 cells 
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remains in the stationary state (Ag, Bo) on the right of the system, so that 
the propagation of the front is not significantly affected by the finite size of 
the system. To draw Fig. 8 (b), we have chosen the parameter values given 
in the caption of Fig. I and imposed k_3 = 8.76 for the following spatial 
profile for the concentration Ro of the inhibitor reservoir: 


Ra = 1, for 1<i< 500, (18) 
Ra = 2.83 x 10741 +0.858, for 500<i< 1000, (19) 
Ra = 1.14, for 1000 <i. (20) 


The simulation is stopped at time teng for which the wave front has passed 
cell i = 1000. It is worth noting that the Turing pattern is unchanged for 
larger values of the final integration time. Then, only the position of the 
concentration gradients associated with the traveling wave evolve in time but 
the Turing structure has stopped growing and remains in a steady state with 
a fixed number of wavelengths. As desired, the increase of the concentration 
Ro leads to the termination of the Turing structure. 

As illustrated in Fig. 7, the Turing structure is expected to be stable 
in the range of 1 < i < 500 for which k_3R = 8.76 and unstable in the 
range of i > 1000 for which k 3Ro = 10. More precisely, according to 
Eq. (13), the maximum of the eigenvalue Umax vanishes for k_3R2 = 9.75, 
i.e. Ro = 1.11 for k_3 = 8.76, which occurs in spatial cell i = 900. Hence, 
the Turing pattern is predicted to be stable in the range of 0 < i < 900 and 
unstable beyond this domain. The results shown in Fig. 8(b) confirm the 
analytical predictions. The amplitude of the spatial oscillations decreases 
between i 2 500 and i ~ 1000. The system is in a steady state in the range 
of 1000 < ż < 1500. 

The increase of Ra not only destabilizes the Turing structure but also 
modifies the steady state values and the propagation speed of the wave front. 
The comparison between Figs. 8 (a) and 8 (b) shows that, as Ro increases, the 
wave front propagates faster, ÅT increases, Br decreases and Bo increases. 
Ås a consequence of the variation of År and Br, the oscillations of Å and 
B concentrations are not symmetrical in the range of 500 < i < 900. The 
decrease of the wavelength predicted in Fig. 7 is more difficult to check by 
a qualitative analysis. Using the numerical results illustrated in Fig. 8 (b), 
we evaluate the local wavelength by computing the number of cells between 
two minima of the Å concentration profile. The results are given in Fig. 9 
and compared to the analytical prediction deduced from Eqs. (11) and (12). 
The agreement between the numerical and analytical results is very satis- 
fying in the range of 600 < i < 900. Oscillations of very small amplitude 
are observed in Fig. 8(b) in the range of 900 < i < 1000, proving that a 
very damped Turing structure remains in a small area where instability was 
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predicted. The wavelength of the structure in the range of 1 < 4 < 500 is 
slightly affected by the increase of Ra from cell i = 500 but the deviation 
from the analytical prediction is only 2.5 percent. This small difference is 
related to the linear approximation used in wavelength evaluation that ne- 
glects non-linear terms that may be more important for large structures. 
Interestingly, the wavelength is sensitively decreased in the expected area in 
which the concentration of the reservoir Ro increases: As shown in Fig. 9, 
the wavelength is reduced from 72 spatial cells to less than 61, before the 
structure disappears. We conclude that an increase in the concentration of 
the reservoir Ro related to the inhibitor B is sufficient to account for the 
destabilization of the Turing pattern associated with a decrease of the wave- 
length. As anticipated by the results given in Fig. 7, according to which an 
increase of Ra decreases the wavelength A and leads to a negative eigenvalue 
Umax around (Ar, Br), we suggest that an appropriate spatial variation of 
Ro can be used in chemical engineering to stabilize the homogeneous steady 
state and induce a termination of the Turing pattern in a growing system. 
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Fig.9. Spatial variation of the scaled wavelength A/Ax of the Turing pattern 
versus cell label x/Ax for the parameter values given in the caption of Fig. 8 (b). 
Symbols: Results deduced from the numerical integration of Eqs. (4) and (5). Solid 
line: Analytical prediction given in Eq. (12). 


4. Conclusion 


In a biomimetic approach, we have addressed the question of the ter- 
mination of a Turing structure in a growing system. A free boundary is 
imposed at the growing part, which ensures that the wavelength of the pat- 
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tern is not perturbed by fixed boundary conditions. After deriving analytical 
expressions for the stability condition and the wavelength of the structure, 
we perform a systematic analysis of the effect of all dynamical parameters 
on the pattern. Apart from the variation of the diffusion coefficient of the 
activator, a well-chosen variation of the dynamical parameters leads to the 
desired behavior, i.e. the simultaneous loss of stability and the decrease 
of the wavelength. In particular, an increase of the effective rate constant 
k_3R2, where k_3 is the rate constant of the reaction injecting the inhibitor 
from the reservoir at the concentration Ro, is associated with a destabiliza- 
tion of the Turing pattern accompanied by a decrease of the wavelength. 

Imposing a spatial variation of the concentration of the reservoir Ro 
turns out to be an appropriate protocol for chemical engineering. However, 
the proposed procedure imposes the total length of the structure but not 
its number of wavelengths. In the framework of developmental biology, for 
example in the case of the growth of the digits or the spine of the verte- 
brates, the termination process has to respect the total number of segments 
for a possible variation in the length of the global structure. Therefore, it is 
necessary to imagine that the system itself is able to count the number of al- 
ready formed segments and to trigger the variation of a parameter leading to 
smaller subsequently formed segments. If the concept of the Turing structure 
is kept in the formation of biological patterns, the presented results could 
be used to suggest such relevant parameters. The local increase of the rate 
constant k_3 that would be activated when a given number of segments has 
already been formed can be straightforwardly proposed. Similarly, the local 
increase of the rate constant ko controlling the autocatalytic production of 
the activator or the local decrease of the rate constant kı or k3, associated 
with the absorption of the activator or the inhibitor by reservoirs, would lead 
to the desired phenomenon. The local decrease of the diffusion coefficient of 
the inhibitor offers an alternative. The nature of the mechanism that would 
trigger such a response of the system when a given number of segments has 
been created remains an open question. 
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111.3 Scaling of Turing patterns 


11.3.1 Summary of the results 


The lack of adaptability to the global size of the system is topical [82, 83, 84, 85, 86, 87, 88] 
and one of the main objections against Turing-based models in morphogenesis. Specifi- 
cally, in somitogenesis, it is expected that vertebra size adapts to the size of the embryo. 
Just before my arrival in the group, B. Nowakowski and A. Lemarchand introduced a 
model, inspired by the Schnakenberg model [89] and the Gray-Scott model [73] capable 
of addressing this issue at the macroscopic scale [36]. They proposed to address the prob- 
lem in the context of molecular crowding, known to lead to non-negligible effects on the 
chemical mechanism [90, 91, 92, 93, 94, 95]. They considered a concentrated system in 
which the variations of the concentration of the solvent cannot be neglected and admitted 


that the reaction scheme presented in Eqs. (III.1-III.3) is modified as follows 


NE JE (111.21) 

2A+B-Ż+3A (111.22) 

B +S == Rg (111.23) 
-3 


They solved the partial differential equations governing the evolution of the concentrations 
and proved that the wavelength of the Turing pattern can be controlled by the deviation 
from the high-dilution limit, roughly defined as the ratio (cą + cp)/cs of the solute 


concentration and the solvent concentration. 


The deviation from the high-dilution limit is more prevalent in smaller systems such 
as biological cells. The challenge I faced was to adapt the model to simulations of particle 
dynamics based on the direct simulation Monte Carlo method presented in Sec. 1.4.2. Even 
if the number of particles in the system varies, one of the constraints of the simulations 
is to keep the total number of simulated particles constant, in order to check the absence 
of bias in the total momentum and kinetic energy. It may imply the creation of ghost 
particles that slow down the simulation. The scheme given in Eqs. (111.21-111.23) does 
not conserve the number of particles and involves reservoir particles that are also time 


consuming. 
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I proposed to solve the problem using the following scheme 


Aas aoe 
ZAJE GA (111.24) 
Bog 496 


Re Sp 


(111.25) 
S —> Rg 


in which the solvent particles play the role of particles of the reservoirs Rg and Rg. In 
other words, when the process given in Eq. (111.25) occurs in a given spatial box, the 
particle B is created with the velocity and the position of a randomly chosen particle S of 
the same box. At the same time, the particle S disappears exactly at the same constant 
rate k_3 as the particle B is created. The step cannot be written S => B, because it 
would introduce a term k 35 in the rate equation of B instead of the constant term k. 3. 
I performed DSMC simulations of a reactive ternary mixture of A, B, and S particles 
with different diameters in order to reproduce different diffusion coefficients as presented 
in Sec. 1.4.2. My results show that the wavelength of Turing patterns can be tuned at 
the submicrometric scale by controlling the total concentration, i.e. the deviation from 
the high-dilution limit. More precisely, doubling the concentration of the solute decreases 
the wavelength of the structure by a factor of 2. The results can be considered as a pos- 
sible interpretation for proportion preservation of embryos in morphogenesis. We suggest 
that they could be used to design biomimetic materials with controlled submicrometric 


properties in chemical engineering. 


11.3.2 Publication 


The results are published in the article “Scaling of submicrometric Turing patterns in 
concentrated growing systems”, G. Morgado, B. Nowakowski, and A. Lemarchand, Phys. 
Rev. E, 98, 032213 (2018) [37]. 
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Chapter IV 


Fisher-Kolmogorov-Petrovsky-Piskunov 


front 


For more than eighty years, the Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) wave 
front has been providing new puzzles to researchers in dynamical systems theory and 
statistical physics. Within the framework of a deterministic description, corrections to 
the asymptotic propagation speed have been determined depending on the steepness of the 
initial condition first by Bramson [96, 97], then by Ebert and Saarloos [98], and currently 
by Brunet and Derrida [99, 100]. The role of fluctuations on the propagation speed has 
been first numerically detected using Langevin equations [101, 102], a master equation, 
DSMC and molecular dynamics simulations [103] The analytical stochastic description of 
FKKP fronts continues to be developed [104]. My contribution to the subject differs by 


taking into account different diffusion coefficients for the two reacting species. 


IV.1 State of the art for identical diffusion coefficients 
The FKPP model generalizes the Verhulst model, presented in Eq. (1.39), 
A+B spå (IV.1) 


to inhomogeneous systems. For species A and B with the same diffusion coefficient D, 


the rate equations are written 


Orca = kcacp + D& ca (IV.2) 
cg = —kcacp + Dcp (IV.3) 


According to Eq. (IV.1), the quantity c*% = cą + cp is conserved, leading to the single 
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cx, cx, 


tt 


> 
w 


Fig. IV.1 Left: Concentrations of A and B species vs. space coordinate z at 
time t = 0. Right: Concentrations of A and B species vs. space coordinate x 
when t oo. The wave travels with a stationary velocity v*. 


rate equation 


c 

Oca = kctcą (i — Sar) + Deca (IV.4) 
c 

This equation exhibits traveling wave solutions, here between the unstable steadv state 


ca = 0 and the stable steady state ca = c. 


IV.1.1 Minimum velocitv v' 


The specific properties of an FKPP front can be qualitativelv understood as follows. The 
unstabilitv of the stationarv state ca — 0 makes the leading edge of the front sensitive 
to perturbations. From a theoretical point of view, the properties of the leading edge in 
which the concentration ca is small can be studied within the framework of a linearized 
analvsis. The front is 'pulled' bv the leading edge and the propagation speed does not 
depend on the nonlinearities of the dynamics. 


The propagation speed v of the front is derived from the linearized version of Eqs. (IV.4) 


fol) 


around the unstable state (ca = 0, cg = c!*) in the moving frame Ç = x — vt 


dca tot deca 
NC = kc” cą + D dc (IV.5) 
Provided that the profile cą follows an exponential form in the leading edge 
cą że % (IV.6) 
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Fig. IV.2 Diagram representing the two branches y,(v) and y_(v) for k = 
tot 
c"=D=1. 


where y depends on v, Eq. (IV.5) leads to a second-order polynomial 
Dy? — vy + ke = 0 (IV.7) 


with real solutions 


VENY v2 = 4kctot D 
L= IV.8 
Y+ 3D (IV.8) 


if v? — 4kcttD > 0. 
For sufficiently steep initial conditions, in particular for the step function shown in 
Fig. IV.1, the wave front converges towards a stationary profile that travels at the mini- 


mum speed v* [105]. According to Eq. (IV.8), the minimum speed v* is given by 
v* = 2V kť* D (IV.9) 


Figure IV.2 shows the two branches 44 (v) and y_(v). The minimum speed v* corresponds 


to the meeting point y+(v*) = y_(v*) of the two branches. 


IV.1.2 Cutoff effect 


In 1997, Brunet and Derrida showed that the introduction of a cutoff in the leading edge 
of the front significantly reduces the minimal velocity of the front. In this section, I recall 
the main lines of the demonstration [106]. If a cutoff e is introduced in the reaction term 


of Eq. (IV.1), the rate equation for concentration c is 


Grea = kctcą (1 — Ada — e) + Dóźca IV.10) 
(tot 


(tot 
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c A 
\c(¢) 
\ 
\ å a 
e(c) ~ 1 | Se <e(t)<1: e(l) <e 
Region \ : Region : Region 


l > $ II r M 


Fig. IV.3 The three regions introduced by Brunet and Derrida [106] in the 
leading edge of an FKPP wave front with a cutoff e. 


where O(x) is the Heaviside step function. Introducing the scale variables 


c= mot (IV.11) 

tot 
pe = ij (IV.12) 
t = ket (IV.13) 


and the coordinate ¢ = 2’ — vt’ in the moving frame, and looking for stationary solutions 
c(C), we find 


ud +c" +c(1— c)O(che)=0 (IV.14) 


Denoting the velocity v- of the front with a cutoff, Brunet and Derrida introduce the shift 


A with respect to the minimum velocity 
A=" = Uz (IV.15) 


As shown in Fig. IV.3, three different regions can be defined in the leading edge: In 
region I, c is of order 1, in region II e < c < 1, and in region III c < e. In region I, the 
front is not significantly affected by the cutoff and the differential equation is expected to 
be the same as Eq. (IV.5). In region II, the concentration c is negligible with respect to 


2 


1 and the nonlinear term c” can be neglected compared to the linear term. In region III, 


the reaction term vanishes due to the Heaviside funtion. 
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It reads 


Regionl v*e+ec’+c=0 
Region II ved +c’ +c=0 (IV.16) 
Region III ved + ec’ = 0 


These three second-order linear differential equations can be solved. The main issue is to 


build a continuous, derivable solution at the boundaries of the different regions. Denoting 


y* and yy + iqi the solutions of Eq. (IV.8) for v* and ve respectively yields 


Region I etc) = Cite 74 
Region II  cyr(6) 2 Cie Y sin (yit + Chr) (IV.17) 
Region III crrr(6) Y ee v-(C-60) 


where Cr, Cir, Cir, and Ġo are constants that can be derived from the boundary condi- 


tions. Between the regions I and II, the boundary condition imposes 
Cite ("7704 = Cyr sin (v6 + Cir) (IV.18) 


On the one hand, it is expected that, according to Eq. (IV.15), the difference y* — y, is 
of order A. On the other hand, Eq. (IV.8) shows that y; is of order AW2. Therefore, 
imposing Cf, = 0 at the leading order in A"? yields 


Cr = Cirvi (IV.19) 


Between regions II and III, the concentration is equal to the cutoff e and ¢ = (9. The 


conditions of continuity and derivability of the function c(¢) are given by 


| Cre~%® sin (Yio) = Eyi (IV.20) 
Cqe 71% [—4, sin (yiCo) + y cos (yido)] = —veeyi 
Combining these two equations gives 

shy W (IV.21) 


— tan(7iCo) 


Intuitively, the difference A is expected to be small, i.e. yy ~ y* = 1 and vz => v* = 2. 


Therefore, it is possible to write 
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which is ensured only if yio ~ 7 + yi. Introducing this last assumption in Eq. (IV.20) 
and assuming that (9 > 1 leads to 


l 
Gx (IV.23) 
Y 
T T Ty* 
¡ x x IV.24 
Bel % a ) 
Brunet and Derrida expand vs into power series of yi, 
m * ga, As so ly * 2 IV.2 
ve = u(y" E iyi) = OY") — gu (PV (IV.25) 
and find that the shift in velocity due to the cutoff obeys 
v! (g*r? 
As —~ IV.26 
2(In e)? ( ) 
py, 
For the same parameter values as in Fig. IV.2, the shift is of order ————. The in- 


2(In e)? 
troduction of a cutoff in the deterministic equation has been shown to correctly repro- 


duce the effect of fluctuations in different stochastic systems that can be associated with 
Eq. (IV.4) in the macroscopic limit. In particular, branching Brownian motion [107] and 
the reaction-diffusion master equation associated with the scheme A + B — 2A [103] 
both lead to corrections of the propagation speed obeying Eq. (IV.26). Qualitatively, 
the discrete nature of the random variables in the two considered stochastic approaches 
implies the existence of a rightmost particle, which plays the role of a cutoff in the leading 


edge of the front. 


IV.2 Results for different diffusion coefficients 


IV.2.1 Deterministic description 
a) High-dilution limit 


The result given in Eq. (IV.9) is obtained for D4 = Dg. I address the more general case 
Da A Dg, which implies that the quantity ca + cg is not constant. It is to be noted that 
the deterministic model with different diffusion coefficients cannot straightforwardly be 
associated with elementary diffusion processes at the particle scale. Indeed, in a binary 
mixture, the diffusion coefficients of the two species are identical as shown in Eq. (1.77). 
A ternary mixture involving a solvent S in addition to A and B species offers a possible 
microsocopic picture of a model with Da 4 Dp., as shown in Sec. 1.3.2. The excess of 


solvent with respect to the solute implies that the collisions between S and A particles, 
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on the one hand, and S and B particles, on the other hand, determine the diffusion 
coefficients of A and B, respectively. Therefore, we introduce a solvent S that allows the 


diffusion coefficients of A and B species to be different. The rate equations are 


ca = kcacp + Daca (IV.27) 
cg = —kcącpg + DpóźcB (IV.28) 


Linearizing the equations in the moving frame Ç = x — vt around (ca = 0,cp = œ) 


leads to 


—VÓĘCA = kea + Daca (IV.29) 
—UOĘCB = —ke'ca + DpóźcB (IV.30) 


where c” is the boundary value of cg on the right side of the system, According to 
Eq. (IV.29), cą does not depend on cg, which means that the same procedure as D4 = Dg 


can be applied. We conclude that the minimum velocity of the front is given by 
u' = 2ykcl DA (IV.31) 


for all Dg. Intuitively, the leading edge asymptotically tends to the state (0, c°), for any 
given value of Dg. The propagation speed of pulled fronts being imposed by the leading 
edge, it is not surprising that the velocity does not depend on Dp. 

In the high-dilution limit, the challenge was to find properties of the front profile 
susceptible to be affected by the difference of diffusion coefficients between species A 
and B. The linearization of the Eqs. (IV.27) and (IV.28) does not help in achieving this 
objective. I worked to develop an analytical approach, important to test the quality of 
the numerical results for possibly small perturbations of front properties. I focused on an 
expansion method proposed by Murray to give an estimation of the profile width of an 
unperturbed FKPP front [2]. The idea implemented by Murray is to consider 1/v? as a 
small parameter. Clearly, the quality of the expansion is not excellent, vx being equal to 
2 for the scaled variables given in Eqs. (IV.11-IV.13). As a consequence, the first-order 
expansion delivered results valid in a small interval of Dp close to the value set for Da 
and I was obliged to determine the second-order corrections. According to the approach 
of Murray, it was legitimate to first consider whether the width of the front would be 
affected by different diffusion coefficients for A and B. Even if the width is perturbed, the 
effect remains small. By examining the results of the numerical integration of Eqs. (IV.27) 
and (IV.28), I then stated that the vertical shift between the profiles of A and B could 
be a good candidate. I proposed to use what I called the height h between the A and B 
profiles, defined as the difference c4(¢ = 0) — cp(ć = 0), where the origin ¢ = 0 of the 
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moving frame is set at c4(¢ = 0) = c*%/2. The results are very satisfying. The height h 


tot 


changes sign for Da = Dg and reaches more than 5% of c’®’ in the investigated range of 


Dg smaller than D4 and more than —25% of ¢* for values of Dp larger than Da. 


b) Concentrated system 


The results I obtained in the dilute case were very encouraging regarding the initial goal of 
using the FKPP front as an indicator of diffusion perturbation in a concentrated system. 


I therefore considered the modified rate equations associated with the FKPP model 


Orca = keace + Dads ( = 4) Oca] — Dpó, ED (IV.32) 
C C 
die = —kcacp + Dpd, (i - c: Ocz — Ds» (Eca) (IV.33) 
C C 


obtained from Eqs. (IV.27) and (IV.28) by taking into account the modified Fick's law 
given in Eqs. (1.61) and (1.62). The same procedure than in Eqs.(IV.29) and (IV.30) is 
applied, but the minimal velocity of the front remains the same as in Eq. (IV.31). My 
initial motivation for studying the propagation of an FKPP front was to exploit its sen- 
sitivity to small perturbations in order to use the front as a sensor of the perturbation 
of diffusion induced by high concentrations. From this point of view, this result, which 
states that the propagation speed does not depend neither on the diffusion coefficient of 
species B nor the concentration of the system, is disappointing. However, the results that 
I obtained within the framework of a stochastic description based on a master equation, 


with the introduction of a cutoff, interestingly challenges the result given in Eq. (IV.31). 


I used the same expansion technique as in the dilute case to determine analytical 
expressions of the width and the height in the concentrated case. The results were con- 
firmed by the numerical solutions of Eqs (IV.32) and (IV.33). The height h proved to 
be a good criterion to reveal the perturbation of diffusion in a concentrated system. The 
high-concentration-induced correction to the height monotonically decreases as Dp in- 
creases and remains larger that 5% for Dg = 20D4. Values of Dg smaller than Da lead 
to high-concentration-induced corrections to the height of more than 25%. However, this 
good score is partly due to the intrinsically small values of h in the range Dg < Da. In 
particular, it should not be forgotten that h vanishes for Dg = Da and cannot be used 
in this specific case. 

To conclude, the experimental determination of the vertical shift h between the profiles 
of the two species at the origin of the moving frame and the comparison with the expected 
value in the high-dilution limit should provide a satisfying test of high-concentration- 


induced perturbation of diffusion, more accurate in the range 2D4 < Dg < 20D4. 
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c) Publication 


The results about the deterministic approach to the perturbation of an FKPP front 
for different diffusion coefficients in the dilute and concentrated cases are published in 
“Fisher-Kolmogorov-Petroskii-Piskunov wave front as a sensor of perturbed diffusion in 
concentrated systems”, G. Morgado, B. Nowakowski, and A. Lemarchand, Phys. Rev. E, 
99, 022205 (2019) [43] 
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The sensitivity to perturbations of the Fisher and Kolmogorov, Petrovskii, Piskunov fron 
quantity revealing perturbations of diffusion in a concentrated solution of two chemical 
diffusivities. The deterministic dynamics includes cross-diffusion terms due to the de n 
limit. The behaviors of the front speed, the shift between the concentration profiles AS species, and the 
width of the reactive zone are investigated, both analytically and numerically. „CJ een the two profiles 


turns out to be a well-adapted criterion presenting noticeable variations with the,dé from the dilution limit 
in a wide range of parameter values. 
DOL: 10.1103/PhysRevE.99.022205 $ © 


O 
© 


© 


2470-0045/2019/99(2)/022205(8) 022205-1 ©2019 American Physical Society 


http://rcin.org.pl 


MORGADO, NOWAKOWSKI, AND LEMARCHAND PHYSICAL REVIEW E 99, 022205 (2019) 


A 
$ 
„dY 
Q 
© 
S 
O 


ER 
> 


Q 


022205-2 


http://rcin.org.pl 


FISHER-KOLMOGOROV-PETROVSKII-PISKUNOV WAVE ... PHYSICAL REVIEW E 99, 022205 (2019) 


022205-3 


http://rcin.org.pl 


MORGADO, NOWAKOWSKI, AND LEMARCHAND PHYSICAL REVIEW E 99, 022205 (2019) 


A 
$ 
„dY 
Q 
© 
S 
O 


ER 
> 


Q 


022205-4 


http://rcin.org.pl 


FISHER-KOLMOGOROV-PETROVSKII-PISKUNOV WAVE ... PHYSICAL REVIEW E 99, 022205 (2019) 


022205-5 


http://rcin.org.pl 


MORGADO, NOWAKOWSKI, AND LEMARCHAND PHYSICAL REVIEW E 99, 022205 (2019) 


A 
$ 
„dY 
Q 
© 
S 
O 


ER 
> 


Q 


022205-6 


http://rcin.org.pl 


FISHER-KOLMOGOROV-PETROVSKII-PISKUNOV WAVE ... PHYSICAL REVIEW E 99, 022205 (2019) 


022205-7 


http://rcin.org.pl 


MORGADO, NOWAKOWSKI, AND LEMARCHAND PHYSICAL REVIEW E 99, 022205 (2019) 


A 
$ 
„dY 
Q 
© 
S 
O 


ER 
> 


Q 


022205-8 


http://rcin.org.pl 


90 Chapter IV. Fisher-Kolmogorov-Petrovsky-Piskunov front 


IV.2.2 Stochastic description 


The results I obtained in Sec. IV.2.1 within the framework of a deterministic description 
show that introducing different diffusion coefficients D4 4 Dg does not modify the wave- 
front speed. 


The contribution of a stochastic description involving discrete random variables, as is 
the case for a master equation, is specially valuable. Using Gillespie algorithm recalled 
in Sec.1.4.1, I performed simulations of the master equation associated with the reaction 
A+B — 2A and Da  Dp. The results are qualitatively different from the results of the 
deterministic description. The master equation predicts that the propagation speed of the 
front is sensitively decreased if species B diffuses faster than species A. Typically, the front 
speed is reduced by 22% for Dg = 16D4. This result can be qualitatively interpreted as 
follows. The fast diffusion of particles B quickly brings them to the vicinity of A particles 
where they are consumed by the autocatalytic reaction. Contrary to intuition, a large 
value of Dg leads to a smoother B profile. Hence, the rightmost particle A is surrounded 
by a smaller number of particles B than in the case Dg < DĄ and the front propagates 
more slowly. 

From a theoretical point of view, this result is the most striking contribution of my PhD 
work. It is rewarding to bring out a new result on a problem that is more than 80 years 
old. 

Using a deterministic analogy inspired from the cutoff approach of Brunet and Derrida, 
I consider Eqs. (IV.27) and (IV.28) in which a cutoff e is introduced 


Oca = keacnO( a — e) + Daca (IV.34) 
vcs = -kcacnO( — e) + DpóźcB (IV.35) 


The linear analysis leads to a correction to front speed obeving Eq. (IV.26) which does 
not account for the behavior at large Dp. The problem is nonlinear and the linear cutoff 
approach presented in Sec. IV.1.2 fails. As suggested in Fig. IV.4, I conjecture that, for 
Dp > Da, the leading edge of the A profile sees a smaller concentration of B than for 
Dp < DĄ, leading to the empirical formula 


T2 
Us = 24 / 2cg.DA ( = xan) (IV.36) 


tot — € is deduced 


where the B concentration cp, at the abscissa z, for which c4(x,)/c 
from the numerical integration of the deterministic equations. 
The results of the master equation satisfactorily agree with the empirical formula 


in which, according to reference [103], the cutoff e is evaluated by the inverse of the 
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number of particles in the reactive interface. Despite intensive efforts and trials involving 
many methods, I could not derive an analytical prediction of the correction to the front 
speed for Dg > Dą. Using the perturbative approach in power of 1/v? developed in the 
deterministic case as explained in Sec. IV.2.1 does not help as the nonlinear effect of Dp 
is lost when Eqs. (IV.34) and (IV.35) are linearized. Deducing the variance (c4cg) from 
a Langevin approach as in Sec. II is of no use because it involves continuous variables 
and even misses the linear cutoff effect. Applying the Hamilton-Jacobi technique to 
solve the master equation [108, 109] is also not successful. The effect of fast diffusion 
of the consumed species on the propagation speed that I numerically evidenced opens 
new perspectives in the fundamental description of Fisher - Kolmogorov, Petrovsky and 


Piskunov wave fronts. 


a) Concentrated system 


In a stochastic description using a master equation, diffusion is a jump process from one 
spatial cell to an adjacent cell. In a dilute system, the diffusion rate of a given species 
only depends on the number of particles of that species in the departure cell. The main 
difficulty in the master equation approach to a concentrated system is the definition of 
the transition rates including cross-diffusion. 

I considered the master equation given in Eq. (1.30) and wrote the associated diffusion 


term as 


P(B), t)|dittusion = PAP ANA — 1) — 1, Na(i) + 1)) 


+ TN i PUN Ali) + 1,N4(i + 1) = 1} 
T Tupo PANB(i m 1) = 1, Ng(i) + 1}) 


+ Tyee (Neli) +1, Na(i + 1) — 1))] (IV.37) 


where TN x (i) are the transition rates associated with the jump of a particle X=A,B from 
cell i containing Nx (i) particles to the left (-) or the right (+), respectively. The transi- 


tion rates must be compatible with the macroscopic diffusion fluxes given in Eqs. (1.61) 


and (1.62). Consequently, TN x (i) has a nontrivial dependence on the particle numbers of 


the two species in both the departure and arrival cells. 


In order to propose an appropriate expression of TN x (y I introduced a discrete flux 
Jx(i + 1/2) at the interface between the cells i and i + 1 and related it to the difference 
of transition rates to the left and to the right. Using Eqs. (1.61) and (1.62) and replacing 

Nx(i + 1) — Nx(i) 
O. 
Bea y NAT 


where (2 stands for the volume of a single cell, I assigned 
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Fig. IV.4 (a) Numbers NĄ of particles A (red dashed line) and Np of parti- 
cles B (black solid line) versus spatial coordinate x deduced from simulation of 
the master equation for Da = 1, Dg = 16 (other parameter values given in 
reference [44]). The vertical dashed line indicates the rightmost cell occupied 
by A particles. (b) Concentrations cą of species A (red dashed line) and cp of 
species B (black solid line) versus spatial coordinate x deduced from numerical 
integration of the deterministic equations in the presence of a cutoff e. The ver- 
tical dashed line indicates the abscissa x, for which the scaled A concentration 
ca(a-)/c’* reaches the cutoff value e. The horizontal line indicates the value cp, 
of B concentration at the abscissa ve. 


well-chosen terms of the flux jy(i + 1/2) to the transition rates to the left and to the right 


Da Na(t + 1/2 : a 
TRAC) Ar? Na(i) ar 2) IDANA(i)— DBNpg(i  1)] (IV.38) 
+ D Np(i + 1/2 : L 
Nel) = 442 NE) vate A 2 [DpNp(i) — DaNa(i + 1)] (IV.39) 


in order to ensure that TH y (3) is positive or equal to zero for any number of particles. 
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The expression of Ti, ,, depends on the number of particles N x(1+1/2) at the interface 
between two cells. Various definitions of this number of particles may be proposed. I 
checked that different definitions, all ensuring that the transition rate vanishes when the 
departure cell is empty, lead to similar results. 

Simulations of the resulting master equation using Gillespie algorithm have been per- 
formed for different values of the diffusion coefficient Dg of species B. I found that the 
decrease of the propagation speed of the front observed as Dp increases is mitigated by 


cross-diffusion which reduces the impact of different diffusion coefficients. 


b) Publication 


In addition to evidencing fluctuation effects on front speed, which constitutes the major 
result, I characterized profile width W and height h in the stochastic approach, in both 
the dilute and concentrated cases. The results about W and h are closer to what could 
be expected after the deterministic study [43]. All my results about the stochastic ap- 
proach to FKPP fronts are published in the article “Stochastic approach to Fisher and 
Kolmogorov, Petrovskii, and Piskunov wave fronts for species with different diffusivities 
in dilute and concentrated solutions”, G. Morgado, B. Nowakowski, and A. Lemarchand, 
Physica A, 558, 124954 (2020) [44]. 
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Conclusion 


In this work, I have been concerned with the formation of structures in living organisms. I 
chose to describe biological systems at an intermediate, mesocopic scale, using methods of 
statistical physics, such as kinetic theory and stochastic methods. My contribution to the 
stochastic elimination of fast variables in the framework of chemical kinetics [29] has illus- 
trated the complex interplay between nonlinear deterministic dynamics and fluctuations. 
Specifically, I have shown that the linearized Langevin equations used to analytically com- 
pute the correlations of concentration fluctuations around steady values do not correctly 
capture the behavior of the system, in particular close to bifurcations. As a perspective, 
with the aim of developing an improved analytical approach, I suggest to consider the 
procedure proposed by Roberts and collaborators [110, 111, 112, 113] to derive stochastic 
normal forms valid in the vicinity of the center manifold. Roberts et al. assign orders of 
magnitude to the different terms of the evolution equations, in particular the Langevin 
forces. They derive approximate equations for the slow variables, describing the dynamics 
on the center manifold up to the desired order. The technique generalizes the ideas devel- 
oped by Arnold to derive deterministic normal forms [114]. The treatment of the vicinity 
of a bifurcation is naturally included in the method and simply consists in including 
small terms related to the bifurcation parameter in the expansion procedure. It is to be 
noted that the expansion leads to non trivial noise terms involving convolution integrals 
of exponential terms and Langevin forces, which introduce memory noise terms in the 
reduced stochastic dynamics. Higher order terms of the expansion even yield nonlinear 
noise combinations between a Langevin force and a convolution integral of an exponential 
function and a Langevin force, which shows how complex the stochastic elimination of fast 
variables is [110]. The quality of the approach, i.e. the order at which the expansion can 
be troncated, could be checked by numerical integration of the Langevin equations, using 
the Euler-Maruyama method with Itó interpretation of the multiplicative noise [115]. The 
great advantage of the numerical integration of the Langevin equations would be to facili- 
tate the comparison with the simulation of the master equation. As proved in the study of 
time asymmetry of correlation functions in far-from-equilibrium conditions in a bistable 


system [63], in oscillating systems close to a Hopf bifurcation [63, 62] and a saddle-node 
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infinite period bifurcation [63], it is however probable that the Langevin approach, al- 
though including multiplicative noise terms and all the nonlinearities of the dynamics, 
will not be sufficient to capture the subtleties of the interaction between fluctuations and 
deterministic nonlinearities. In general, the non Gaussian character of chemical noise 
intrinsically invalidates the Langevin approach. It is especially true in the intricate case 
of the stochastic elimination of fast variables. The correct description of concentration 
fluctuations is essential to predict the behavior of many systems that are sensitive to 
perturbations, in particular due to explosive behavior, the vicinity of a bifurcation, the 
existence of many simultaneously stable states or simply because they are small [116]. 
Hence, combustion hazards, pattern formation in developmental biology including Turing 
structures and Fisher - Kolmogorov, Petrovsky, Piskunov (FKPP) wave fronts require 
stochastic analyses. 

The simulations of a submicrometric Turing pattern in a concentrated system I performed 
refute certain objections to Turing’s model regarding the preservation of proportions in 
embryos. Assuming an appropriate role of the solvent in the chemical mechanism is suf- 
ficient to control the wavelength of the structure by monitoring the concentration of the 
solution. A significant decrease of the wavelength is obtained in a more concentrated so- 
lution for the considered chemical scheme. The adaptation of the size of the structure to 
the size of the embryo then follows from the hypothesis that a smaller embryo has greater 
concentrations of morphogens, imposed by the mother and not by the size of the embryo. 
In this respect, small embryos lead to a crowded environment [117, 118, 119]. The re- 
sults can be exploited to design materials with controlled submicrometric properties in 
chemical engineering [120, 121, 122, 123]. Following a biomimetic approach, I proposed 
experimental conditions, compatible with the requirements of chemical engineering, to 
observe the termination of a Turing structure in a growing system. Among the different 
parameters playing on the stability of the pattern and the value of the wavelength, I 
selected the concentration of the reservoir which sets the injection rate of the inhibitor 
into the system for its easy control in the region where the experimentalist wishes the 


structure to stop. 


I also characterized the effect of concentrated solutions on another spatio-temporal 
structure, often encountered in biology, a Fisher - Kolmogorov, Petrovsky, Piskunov 
(FKPP) wave front, used to model the propagation of a virus or a favorable genetic trait 
in a population. I focused on the modification of diffusion due to high concentrations. 
Beyond the description of concentrated solutions, I was led to consider reactive species 
with different diffusion coefficients. Performing simulations of the master equation, I ob- 
tained a nice result stating that an FKPP front is slowed down when the consumed species 
diffuses faster than the autocatalytically produced species. The analytical description of 


this phenomenon, using a nonlinear cutoff approach, is certainly a reasonable perspec- 
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tive. The consequences of high concentrations on biological structures deserves further 
attention. The results show that, using particle dynamics simulations, it is possible to 
check the validity of macroscopic models at the submicrometric scale. The new directions 
explored in the field of chemical engineering, as part of a biomimetic approach, could 
encourage experimentalists to design model systems for testing the main influence of high 


concentrations on organized behaviors. 


The results obtained during my PhD led me to formulate some general ideas about the 
modeling of structures in biology and prompted me to revisit how the vision of biological 
patterns had evolved over the last century and the beginning of 21st century. Morphogen- 
esis is an important part of developmental biology. Axial segmentation and the formation 
of periodic patterns are observed in invertebrates such as insects and crustaceans as well 
as in vertebrate embryos. The striking analogies between biological structures and the 
patterns that spontaneously emerge in far-from-equilibrium chemical systems logically 
incited theoreticians to use models of chemical kinetics to study biological phenomena. 
Experimental evidence of periodic spatial structures in a reaction-diffusion system, with- 
out the contribution of physical phenomena such as gravity and mechanical instabilities, 
was given in 1990 [72], long after the corresponding model was proposed by Turing in 1952. 
Previously, it had taken even more time to apply the concepts of dynamical systems, de- 
veloped by mathematicians such as Henri Poincaré around 1900, to the description of 
far-from-equilibrium nonlinear phenomena in uniform systems. The temporal organi- 
zation in homogeneous chemical systems, from periodic oscillations to chaos, has been 
interpreted in the context of dynamical systems theory in the 70s, first from a theoretical 
point of view, in particular in the group of the Nobel prize Ilya Prigogine in Brussels. 
Then, in the 80s, the Belousov-Zhabotinsky (BZ) reaction has been extensively stud- 
ied. With the development and mastery of continuously stirred tank reactors (CSTRs) 
ensuring far-from-equilibrium conditions, the BZ reaction proved to provide an ideal ex- 
perimental example of dynamical systems exhibiting all the different types of bifurcations 
and scenarios to chaos [124, 125, 126, 127, 128]. Application to biology arrived later and, 
at the very beginning of the 2000s, the concept of systems biology was introduced. It is 
indeed tempting to apply the notions of the theory of dynamical systems to model cer- 
tain biological phenomena. Living systems are typically maintained far from equilibrium 
by exchanges of energy and matter with the environment and spontaneously evolve into 


organized structures. 


Two antagonistic views, known as holism and reductionism, are usually adopted to 
study chemical systems and provide disjointed information. On the one hand, holism 
proposes a global approach to a system and neglects the discrete nature of matter. After 
a modeling effort to identify essential mechanisms and extract a small number of dynam- 


ical variables, the system is described at the macroscopic scale. It is then possible to 
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write differential equations governing the evolution of a homogeneous system and partial 
differential equations in the case of an inhomogeneous system. The nonlinearities of the 
equations determine the behavior of the system, depending on the parameter values and 
the involved bifurcations. A global, general classification of some biological phenomena 
can be obtained within the framework of such a deterministic, macroscopic analysis. Such 
a description ignores the fluctuations of the macroscopic variables that are induced by 
the microscopic dynamics of a huge number of elementary constituents. This global ap- 
proach can be misleading in systems of small size like the systems typically studied in 
developmental biology. The role of fluctuations on the macroscopic behavior of a system 
is also significant, even in a large system, close to bifurcations [29] or for marginally stable 
solutions such as FKPP fronts [44]. 


On the other hand, reductionism focuses on the atomistic or molecular scale. However, 
the huge number of particles prompts the use of numerical simulations. The evolution of 
wave functions or particle positions and velocities, dictated by fundamental interactions 
according to quantum or classical mechanics, can then be followed. Even with the lat- 
est generation of computers, it is still extremely long to reach the space and time scales 
necessary to observe the formation of macroscopic patterns using ab-initio simulations, 
density functional theory, and even molecular dynamics. Moreover, the connection be- 
tween changing a parameter characterizing atomistic interactions and changing a macro- 
scopic property of a structure is not firmly established. The numerical simulations at 
the microscopic scale may provide empirical knowledge but it remains difficult to build 
qualitative relationships between the molecular structure and the macroscopic properties 
of the system. Reductionism gives access to the specific properties of a given system but 
the prediction of the behavior of a system showing a small difference in the microscopic 


characteristics often fails. Generalization within a reductionist approach is delicate. 


In order to bridge the gap between the two remote microscopic and macroscopic scales, 
I chose to develop a description at the mesoscopic scale, which presents the advantage 
of accounting for the fluctuations as well as providing a general framework to develop 
analytical approaches. The strength of statistical physics is to propose a probabilistic de- 
scription of a system composed of a large number of components. Microscopic features of 
molecules such as their atomic components, their chemical functions, and their structure 
are not explicitly taken into account by the stochastic description, which retains some 
consequences of these microscopic properties at the mesoscopic scale and take them into 
account in a coarse-grained manner. In particular, the notion of stochastic trajectory 
includes dissipation and irreversibility but with a more refined approach than the macro- 
scopic description. Specifically, microscopic reversibility with respect to time reversal is 
lost in a master equation approach [61, 62, 63] but the production of entropy may de- 


crease along a stochastic trajectory. The fluctuation theorem demonstrated by Galavotti 
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and Cohen makes precise the probability that such an event occurs [129, 130, 131, 132]. 
Writing a master equation associated with a given system, for example to account for 
cross-diffusion in a concentrated system [44], already represents an effort of modeling, i.e. 
an effort to extract ingredients from the microscopic scale that are essential to describe 
the phenomenon of interest at the mesoscopic scale. Approximating a master equation by 
Langevin equations makes it possible to derive analytical expressions that better point out 
the discontinuities in the properties induced by bifurcations [29]. From a general point of 
view, the analytical approaches made available by statistical physics offer the possibility 
to better « understand » complex phenomena. They at least guide the intuition and make 
it possible to guess the behavior of an unknown system having some points in common 


with the one previously studied. 


The categorization of a given description into a holistic approach or a reductionist 
approach is not always clear. Specifically, the stochastic description using Langevin equa- 
tions is obtained by adding a noise term to the deterministic equations, which includes a 
part of the microscopic complexity into the equations. In this respect, the description by 


Langevin equations lies between the two approaches. 


Kinetic theory offers another interesting view by considering both deterministic macro- 
scopic behavior and intrinsic probabilistic behavior of matter: It looks at the evolution 
of the distributions of the position and velocity of the particles. In particular, collision 
integrals govern the evolution of particle velocities. I have illustrated how kinetic theory 
is able to make the link between the elastic and reactive collisions between hard spheres 
and macroscopic parameters such as diffusion coefficients and rate constants [37]. The 
need for considering a mixture of at least three different species to obtain different diffu- 
sion coefficients for two reactants becomes obvious in kinetic theory whereas it is hidden 
in the two partial differential equations governing the macroscopic evolution of the con- 
centrations of the two reactants [37]. As a conclusion, I would like to point out the value 
of describing a biological system on a mesoscopic scale in order to extract minimal in- 
gredients, decipher mechanisms, and obtain some analytical results guiding intuition and 
understanding [133, 134, 135, 136]. However, the master equation and the kinetic equa- 
tions are not usually solvable without strong approximations. It is of primer importance 
to check the validity of the analytical results. In this respect, both approaches benefit 
from the possibility to perform direct simulations of the equations using kinetic Monte 
Carlo simulations. Gillespie algorithm simulates stochastic trajectories obeying the mas- 
ter equation and the direct simulation Monte Carlo method (DSMC) introduced by Bird is 
a direct simulation of the kinetic equations including fluctuations. In a DSMC simulation, 
a particle represents thousands of true molecules [18], is assimilated to a hard sphere, and 
its collisions are randomly treated in an approximate way. Under these assumptions, par- 


ticle dynamics is simulated up to thousand times faster than using molecular dynamics, 


http://rcin.org.pl 


110 Conclusion 


which allows us to simulate Turing patterns in a reasonable computation time, even in a 
concentrated system [37]. Finally, it is worthwhile to notice that the efficient simulation 
tools provided by kinetic Monte Carlo algorithms, that I first applied to materials sci- 
ence, specifically to the submicrometric simulation of gypsum crystallization [137, 138], 
are particularly well adapted to the simulations of biological submicrometric systems, 


which gives an idea of the potential of such methods. 
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